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ABSTRACT 


A  documentation  of  a  six-degree-of-freedom  model  for  digital 
simulation  of  th'*  trajectory  of  an  unguided,  fin-stabilized,  wind- 
sensitive  rocket  is  presented.  This  model  was  developed  by  the  At¬ 
mospheric  Sciences  Laboratory,  White  Sands  Missile  Range,  New  Mexico, 
to  study  both  theoretical  and  empirical  performance  characteristics 
of  unguided  rockets. 

The  basic  equations  of  motion  and  their  mathematical  formula¬ 
tion  for  this  model  are  presented  without  derivation. 

A  general  flow  chart,  a  listing  of  the  program,  a  list  of  the 
principal  flads  used,  and  a  listing  of  a  typical  input  data  deck 
are  incited. 
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INTRODUCTION 


In  recent  years  many  theories  anti  computer  programs  have  been 
developed  to  simulate  missile  trajectories.  These  models  range  from 
the  extremely  simple  to  the  very  complex,  the  degree  of  complexity 
usually  depending  upon  the  specific  simulation  requirement  placed  upon 
the  developer. 

This  report  describes  a  six-dcgree-of-frcedom  computer  program 
developed  by  the  Atmospheric  Sciences  Laboratory,  White  Sands  Missile 
Range,  New  Mexico,  to  simulate  the  trajectory  of  an  unguided,  fin- 
stabilized,  wind-sensitive  rochet.  The  simulation  theory  upon  which, 
the  program  is  based  is  presented  in  reference  1.  The  references 
may  be  consulted  for  derivations  of  the  equations  and  coordinate  trans 
formations;  however,  the  coordinate  systems  used  and  the  principal 
equations  which  are  evaluated  by  the  program  are  included  for  complete 
ness. 


The  system  is  programmed  in  Fortran  IV  language.  It  consists 
of  a  main  program  (monitor),  two  subprograms  which  serve  as  submon¬ 
itors  for  specific  simulation  options,  and  a  group  of  subroutines, 
each  designed  for  a  specific  task. 

An  attempt  has  been  made  in  this  report  to  present  a  reasonably 
complete  program  documentation  without  boring  the  reader  with  trivia. 
With  this  goal  in  mind  most  of  the  routines  are  documented  in  four 
parts:  (1)  a  statement  of  the  purpose  of  the  routine  and  the  equa¬ 
tions  to  be  evaluated;  (2)  definition  of  the  principal  flads;  (3) 
a  macro  subroutine;  and  (4)  a  listing  of  the  program  instructions. 

It  will  be  observed  that  in  some  of  the  minor  routines  (2)  and/or 
(3)  have  been  deleted. 


COORDINATE  SYSTEIS  AND  TRANSFORMATIONS 


Three  right-hand  coordinate  systems  are  used  in  the  program. 
These  are:  (1)  The  launcher  system  (denoted  X',  Y',  Z')  which  has 
its  origin  at  the  launcher  and  rotates  with  the  earth.  The  positive 
X'-axis  points  east;  the  positive  Y'-axis  points  north;  and  the  pos¬ 
itive  Z'-axis  along  the  outward  normal  of  the  earth.  The  X'-Y'  plane 
is  tangent  to  the  earth  at  the  launcher.  (2)  The  inertial  system 
(X,  Y,  Z)  has  its  origin  at  the  center  of  the  earth.  The  system  is 
oriented  so  that  the  X-Y  plane  lies  in  the  earth's  equatorial  plane 
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with  the  positive  Y-axis  pointing  initially  through  the  longitude  of 
the  launcher.  The  Z-awis  is  coincident  with  the  earth's  axis  of  ro¬ 
tation  and  is  positive  toward  the  North  Pole.  This  system  uoes  not 
rotate  with  the  earth.  (3)  The  body  system  (x,  y,  z)  has  its  origin 
at  the  center  of  gravity  of  the  rocket.  The  x-axis  coincides  with  the 
longitudinal  axis  of  the  rocket  and  is  positive  toward  the  nose. 

The  position  of  the  y-  and  z-axes  is  determined  by  the  rocket's  mo¬ 
tion.  The  initial  positions  of  these  axes  are  defined  as  follows. 

Let  0  be  the  angle  between  the  x-axis  and  the  positive  Z'-axis  meas¬ 
ured  from  the  Z'-axis.  The  y-axis  lies  in  the  X'-Y'  plane  and  is 
positive  in  the  direction  of  positive  5.  It  is  easy  to  see  that  the 
angle  between  the  positive  Y'-  and  y-axes  is  the  launch  azimuth  meas¬ 
ured  clockwise  from  the  Y'-axis. 

The  transformation  between  any  pair  of  these  coordinate  systems 
depends  upon,  among  other  things,  the  assumed  earth  model.  The  pro¬ 
gram  considers  the  earth  to  be  an  oblate  spheroid  with  an  equatorial 
radius  of  20,926,428  feet  and  an  eccentricity  of  .00672267.  The  trans 
formation  from  the  launcher  system  to  the  inertial  system  is  easy  to 
obtain  from  the  geometry  of  the  problem.  This  transformation  is,  in 
matrix  form, 


— - \ 

~<  X 

j 

* 

L2J 

-  cos  wt 

-  sin  wt 

0 


sin 

A 

g 

sin 

wt 

-  cos 

A 

g 

sin  wt 

) 

V 

sin 

A 

g 

cos 

wt 

cos 

A 

g 

cos  wt 

1 

Y' 

cos 

A 

sin 

A 

Z' 

6  S 


where  A  is  the  geodetic  latitude  of  the  launcher,  w  is  the  earth's 
rotatioS  speed  and  t  is  time  after  launch.  The  derivation  of  the 
preceding  transformation,  except  for  the  assumption  of  a  spherical 
earth,  is  presented  in  reference  1. 

The  transformation  from  the  body  system  to  the  inertial  system 
is  obtained  by  integrating  the  derivatives  of  the  elements  of  the 
transformation  matrix  --  the  direction  cosines.  This  transformation 
is  denoted  by 
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It  it  show-,  in  [1]  that  the  required  derivatives  are 


where  p,  q,  r  are  the  x,  y,  z  components  of  the  rotation  of  the  body 
system  with  respect  to  the  inertial  system. 

It  is  shown  in  reference  1  that  the  initial  conditions  for  these 
direction  cosines  arc 
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where  0  is  the  launcher  elevation  angle  and  a  is  the  launcher  azimuth 
o  *  o 

angie  measured  clockwise  from  north. 


Tllli  MAIN  PROGRAM  LRBM 


LRBM  is  used  to  monitor  the  other  routines,  to  initialize  the 
trajectory  simulation,  and  to  monitor  the  type  of  trajectory  desired. 


INITIALIZATION 


Initially  the  x-axis  and  Z'-axis  form  an 
elevation  of  the  x-axis  angle  and  is  measured 
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the  projection  of  the  x-axis  in  the  X'Y'  plane  forms  an  angle  aQ  with 

the  Y'-axis  (this  is  the  azimuth  angle.,  of  the  x-axis  and  is  measured 
clockwise  from  the  Y'-axis).  Since  the  y-axis  lies  initially  in  the 
X'Y'  plane,  it  forms  an  angle  aQ  +  90  with  the  Y'-axis.  (See  Figure 
1).  Thus  initially,  c 


sine  cosa 
o  o 

-sina 

o 

cosa  cose 
o  o 


X  «  -  X' 

Y  ■  -  Y*  sinX  +  Z'  cosX 

g  g 

Z  e  -  Y'  cosX  ♦  Z*  sinX 

S  g 

X  =  i.  U  ♦  wX 
1  o 

Y  ■  U  ♦  wY 

i  o 

Z  •  l.Un 

3  O 


cose 

o 

0 

-sine 

o 


X' 

Y* 

Z* 


where  X^  is  the  geodetic  latitude  of  the  launcher,  Uq  is  initial  x 

component  of  velocity  in  body  system,  and  w  is  the  angular  rate  of 
the  earth's  rotation. 

The  initial  conditions  for  nu,  ik,  i  ■  1,  2,  3,  are  given 
by 
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NTYPE  is  an  indicator,  road  in  the  input  routine,  that  determines 
the  type  of  trajectory  to  be  simulated  (1,  a  regular  trajectory;  2, 
an  anguiar  displacement  table;  3,  iterative  type  trajectory;  4,  the 
ballistic  factor  option;  or  S,  a  parameter  variability  type  trajec¬ 
tory).  LRBM  us<-.,  NTYPE  to  monitor  the  program  and  obtains  the  desired 
type  of  trajectory  simulation. 

The  equations  of  motion  of  an  unguided  rocket  are  numerically 
integrated  by  a  Runge-Kutta  integration  technique.  This  technique 
is  discussed  in  the  section,  the  Integration  Routines.  LRBM  sets 
up  the  initial  integration  interval  for  each  phase  and  monitors  the 
integration  henceforth. 

Several  arrays  used  by  all  routines  are  defined.  They  are  the 
Y,  DY?  and  R0FF1  arrays.  All  three  arrays  have  double  subscripts. 

The  second  subscript,  J,  denotes  the  step  in  the  Runge-Kutta  integra¬ 
tion.  The  Y  array  then  has  the  following  definitions: 

Y(1,J)  -  time  of  the  trajectory  simulation. 

Y(2,J)  -  x  component  of  velocity  vector  in  body  s>-.ten. 

Y(3,J)  -  y  component  of  velocity  vector  in  body  system. 

Y (4, J)  -  z  component  f  ve.’  city  vector  in  body  system. 

Y (5, J)  -  x  component  of  the  rotation  of  body  system  with  re¬ 

spect  to  the  inertiai  system. 

Y(6,J)  -  y  component  of  the  rotation  of  the  body  system  with 
respect  to  the  inertial  system. 

Y(7,J)  -  z  component  of  the  rotation  of  he  body  system  with 
respect  to  the  inertial  system. 

Y(8,J),  Y (9, J) ,  Y (10, J)  -  direction  cosines  of  the  body  X-axis 

with  respect  to  the  inertial  axes. 

Y(il,J),  Y(12,J),  Y(1 u, J)  -  direction  cosines  of  the  body 

Y-axis  with  respect  to  the  iner¬ 
tial  axes. 

Y(14,J),  Y(15,J),  Y(16,J)  -  direction  cosinta  of  the  body 

Z-axis  with  respect  to  the  iner¬ 
tial  axes. 


Y(i7,J)  -  X  component  of  the  inertial  position. 

Y(18,J)  -  Y  component  of  the  inertial  position. 

Y(19,J)  -  2  component  of  the  inertial  position. 

Each  DY(I,J),  I  equals  2  through  19,  is  the  derivative  of  the 
corresponding  Y(I,J)  entry. 

Each  ROFFI(I,J),  I  equals  2  through  19,  is  the  round-off  error 
associated  with  the  corresponding  Y(I,J)  entry  due  to  the  Runge-Kutta 
integration.  The  principal  flads  are: 

AL  -  Azimuth  angle  of  rocket  on  launcher. 

Bill  -  Current  height  for  cutoff  of  wind  profile  for  ballistic 

factor  option. 

CLAT  -  Cosine  of  geocentric  latitude  of  launcher. 

CLATG  -  Cosine  of  geodetic  latitude  of  launcher. 

CNTR  -  Number  of  the  previous  trajectory. 

DELC  -  Cross  unit  wind  effect  (m/mph) 

DELH  -  Head  unit  wind  effect  (m/mph) 

DELT  -  Tail  unit  wind  effect  (m/mph) 

DONE  -  Indicates  if  integration  should  be  continued  (1)  or 

not  (2). 

DYSTOR(22)  -  Storage  of  DY  array  for  future  use. 

ENF  -  Indicates  whether  the  winds  in  the  subroutine,  AOPTUN, 

have  been  iterated  for  (1)  or  not  (0), 

t 

EPSQ  -  Square  of  the  eccentricity  of  the  earth. 

GCONI,  GNU  -  Constants  used  in  evaluation  of  gravity  for  geodetic 
earth. 

H  -  Integration  interval  (sec). 
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-  Heights  in  wind  tabic  (ft). 

-  Last  entry  used  in  atrr  ,pheric,  wind,  Mach,  and  tine 
tables. 

-  Indicates  if  all  wind  profiles  of  ballistic  factor  op¬ 
tion  have  been  sinulated  (1)  or  not  (0). 

-  Indicates  print  out  only  at  end  of  each  phase  (1). 

-  Indicates  whether  to  call  ballistic  factor  option  (2) 
or  not  (1). 

-  Indicates  print  out  (2)  or  not  (1),  end  of  phase  (3), 
or  impact  (4). 

-  Entry  of  wind  currently  being  used  in  ballistic  factor 
option. 

-  X  and  Y  entry  in  table  for  angular  displacement  table. 

-  Phase  currently  being  used. 

-  Phase  at  which  to  pick  up  trajectory  for  ballistic  factor 
option. 

-  Phase  end  of  which  booster  drops  off. 

-  Total  number  of  phases  for  trajectory. 

-  Indicates  whether  integration  is  accepted  (1)  or  not 

(2). 

-  Indicates  whether  a  regular  trajectory  (1),  iterative 
type  trajectory  (2),  an  angular  displacement  table 
(3),  the  ballistic  factor  option  (4),  or  a  parameter 
variability  trajectory  (5)  is  to  be  simulated. 

-  Indicates  if  the  height  of  the  simulated  trajectory 
is  less  than  or  equal  to  height  of  impact  area  (1) 
or  not  (0) . 

-  Print  interval  (sec). 
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PSL 

RANGE 
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-  Atmospheric  pressure  at  sea  level  (lbs). 

-  Distance  to  desired  impact  point,  (ft). 

-  Radius  of  earth  at  equator  (ft). 

-  Magnitude  of  radius  vector  of  oblate  earth  along  ra¬ 
dius  vector  to  rocket. 

-  Sum  of  RO  and  height,  above  sea  level  of  impact  area. 

-  Storage  of  ROFF1  array  for  future  use. 

-  Sine  of  geocentric  latitude  of  launcher. 

-  Sine  of  geodetic  latitude  of  launcher. 

-  Time,  end  of  current  phase. 

-  Tilt  angle  of  rocket  on  launcher. 

-  Current  time  of  simulated  trajectory  (sec). 

-  Time  rocket  leaves  launcher  (sec). 

-  Time  of  end  of  last  phase. 

-  Tower  tilt  effect  (ft/radian). 

-  Initial  body  coordinates  of  aerodynamic  velocity. 

-  Angular  rate  of  rotation  of  earth. 

-  Theoretical  range  displacement  due  to  wind  (ft). 

-  E-W  and  N-S  components  of  wind  in  launcher  coordinate 
system  (mph) . 

-  Theoretical  E-W  and  N-S  components  of  wind  displace¬ 
ment  of  rocket  (ft). 

-  egocentric  latitude  of  launcher. 

-  Geodetic  latitude  of  launcher. 

-  Total  number  of  trajectories. 


SJSCALL  Y 
AND  DY 
AtiAYS 


XUW,  YUW  -  Cress  and  range  Unit  wind  effects  (m/mph) . 

XWANT,  YWANT  -■  E-lv  and  N-S  components  of  desired  impact  point  (ft). 
ZL  -  Height  f  impact  area  above  sea  level  (ft). 


THE  TWO  SUBMONITORS 


For  convenience  of  the  user,  techniques  are  available  for  simu¬ 
lating  a  series  of  trajectories  from  one  computer  submission.  The 
first  of  these  techniques  simply  consists  of  submitting  a  number  of 
consecutive  data  decks  -  one  for  each  trajectory.  Two  other  techni¬ 
ques  are  available  for  computing  a  series  of  trajectories  of  a  spe¬ 
cialized  nature.  These  computations  are  controlled  by  the  submoni¬ 
tors,  BALFAC  and  AOPTUN. 

SUBROUTINE  BALFAC  (INTA) 

This  subroutine  is  used  to  monitor  the  calculations  of  ballistic 
factors  (wind  weighting  factors)  and  unit  wind  effects.  A  complete 
trajectory  is  simulated  for  each  entry  in  the  wind  table  (see  section 
on  Input).  Each  card  in  the  wind  table  has  an  entry  of  wind  and  height. 
The  trajectory  corresponding  to  a  given  entry  in  the  wind  table  is 
computed  for  a  wind  profile  having  the  wind  value  of  the  entry  when 
the  simulated  altitude  is  below  the  specified  height  and  a  value  of 
zero  when  the  simulated  altitude  is  above  this  height.  The  first  en¬ 
try  in  the  wind  table  should  have  zero  winds,  i.e.,a  no-wind  trajec¬ 
tory  is  required. 

The  X  and  Y  components  of  impact  for  each  entry  are  subtracted 
from  the  no-wir.d  impact  point  to  obtain  a  range  displacement.  The 
range  displacement  of  each  trajectory  is  then  divided  by  the  range 
displacement  of  the  final  trajectory  to  obtain  a  ratio  or  percentage 
of  the  total  range  displacement.  Each  ratio  is  then  subtracted  from 
the  previous  ratio  to  give  the  ballistic  factor  for  that  layer. 

The  unit  wind  effect  is  obtained  by  dividing  the  last  range  dis¬ 
placement  by  the  magnitude  of  the  wind. 

Three  indicators  are  used  by  the  subroutine  to  control  the  logic 
flow.  These  are  INTA,  IKIND  and  I  SPEND. 
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Tlie  routine  is  entered  after  the  completion  of  each  integration 
step.  INTA  determines  the  entry  point  into  the  subroutine.  Muring 
the  simulation  of  a  given  trajectory  the  value  of  INTA  is  1.  'I'.licn 
impact  is  detected,  by  LRU'!,  INTA  is  set  equal  to  2  and  LALFAC  is 
entered  to  store  the  impact  point  and  establish  the  initial  conditions 
for  the  next  simulation. 

I KIND  is  the  indicator  to  determine  which  wind  table  entry  to  use 
for  a  given  trajectory. 

IliPl.ND  is  an  indicator  which  tells  the  routine  if  the  simulations 
have  been  completed.  This  determination  is  made  in  LRU'1.  If  the 
simulations  have  been  completed,  IBPLND  is  set  equal  to  1  and  UALPAC 
is  entered  to  c  aputc  t lie  ballistic  factors  and  unit  wind  effect. 

SUUROIJTINL  AOPTUN  (K) 

The  simulation  model  contains  a  technique  for  determining  the 
launcher  settings  necessary  to  compensate  for  a  given  wind  profile 
so  that  the  computed  impact  will  be  in  a  desired  area.  The  launcher 
settings  arc  determined  by  an  iterative  procedure  which  is  discussed 
in  reference  3. 

This  subroutine  is  used  to  monitor  the  iteration  procedure.  In 
addition  it  is  used  to  monitor  the  computation  of  a  nomogram  of  launcher 
angles  versus  ballistic  winds.  (See  reference  These  latter  cal¬ 
culations  arc  referred  to  as  angular  displacement  trajectories. 

£  is  an  indicator. used  to  determine  the  entry  point  to  the  sub¬ 
routine  and  has  admissible  values  of  1  and  3.  When  K  is  equal  to 
1  the  initial  approximation  of  the  launcher  setting  is  made  and  the 
wind  profile  is  set  up.  The  subroutine  is  entered  witli  a  value  of 
K  =  3  at  impact  and  a  check  is  made  to  determine  if  the  simulated 
impact  is  within  a  required  tolerance  of  the  desired  impact.  If  it 
is  not,  a  new  set  of  launcher  angles  is  determined  am!  the  iteration 
is  continued.  If  the  tolerance  is  met  and  an  iterative  trajectory 
is  being  computed,  a  return  is  generated.  If  a  nomogram  is  being 
generated,  the  launcher  setting  is  stored  in  tabular  form  and  the 
next  wind  profile  is  set  up  unless  the  nomogram  has  been  completed; 
whereupon  a  return  is  generated.  The  first  approximation  for  the 
launcher  settings  for  this  profile  is  made  and  the  new  trajectory 
is  begun.  The  principal  flads  for  the  two  submonitors  arc: 

AI.  -  Azimuth  angle  of  rocket. 

ALPHA  (11,11)-  Storage  of  azimuth  angles  in  tabular  form  for  angular 
displacement  option. 
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-  Ballistic  factors  of  the  rocket. 

-  Ballistic  factor. 

-  Height  of  current  entry  in  wind  tabic  (feet). 

-  Cross  unit  wind  effect  (m/mph). 

-  Head  unit  wind  effect  (m/mph). 

-  Ballistic  wind  displacement  of  current  trajectory 
(miles) . 

-  Tail  unit  wind  effect  (m/nph) . 

-  C-l;  ballistic  wind  displacement  of  current  trajectory 
(miles) . 

-  N-S  ballistic  wind  displacement  of  current  trajectory 
(miles) . 

-  Indicator,  set  equal  to  one  when  angular  displacement 
table  finished. 

-  Heights  used  in  the  wind  table  (feet). 

-  Height  of  rocket  above  oblate  earth  (feet). 

-  Indicates  whether  or  not  the  entire  wind  table  has 
been  used. 

-  Indicates  whether  to  use  wind  or  not. 

-  Entry  in  wind  table  currently  being  used. 

-  Number  of  entries  in  ballistic  factor  table. 

-  Indicates  type  of  trajectory  to  be  simulated. 

-  Number  of  entries  in  wind  table. 

-  Same  ns  DFLR. 

-  Percentage  of  wind  effect. 
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SET  UP  INITIAL  OBTAIN  BALLISTIC 

WINDS 

WIND  PROFILE  FOR  ITERATION 


IS  NTYPE;2? 


STORE  LAUNCHER 
SETTING  IN 
TABULAR  FORM 


IS  TABLE 
COMPLETE 


e  IS  IMPACT  > 
WITHIN  TOLERANCE, 


ITERATE  FOR 

LAUNCHER 

SETTINGS 


RETURN 


RETURN 


SET  UP  NEXT 
WIND  PROFILE 


DETERMINE 

THEORETICAL 

WIND 

DISPLACEMENT 


IS 


»  4*;sgJfi(PKis '  m&jm*0ir* *m*t!*u&#*  >  >■  ■  >•  .  „ 

> 


SUMX  -  E-W  component  of  ballistic  wind. 

SUMY  -  N-S  component  of  ballistic  wind. 

TESTVA  -  Range  difference  between  simulated  and  desired 

impact. 

RH  -  Tilt  angle  of  rocket. 

THETA  (11,11)  -  Storage  of  tilt  angles  in  tabular  form  for  angular 

displacement  option. 

TOWTIL  -  Tower  tilt  effect  (ft/rad). 

UNIT  -  Unit  wind  effect  (m/mph) . 

WIND  -  Magnitude  of  wind  (mph) . 

KXr(50),  WYP(SO)  -  Entries  of  E-W  and  N-S  components  of  wind. 

-  Value  of  E-K  and  N-S  components  of  wind. 

-  E-W  and  N'-S  components  of  difference  between  sim¬ 
ulated  and  desired  impact. 

-  Increments  of  E-W  and  N-S  component  of  wind  for 
angular  displacement  option. 

-  Initial  and  terminal  values  of  E-W  component  of 
wind  for  angular  displacement  option. 

-  Cross  and  range  component  of  unit  wind  effect. 

-  E-W  and  N-S  components  of  desired  impact. 

-  E-W  and  N-S  components  of  simulated  impact. 

THE  EVAL  ROUTINE 

The  subroutine  EVAL  is  used  to  evaluate  the  rocket  equations 
of  motion  and  to  compute  the  derivatives  of  the  direction  cosine  ma¬ 
trix  which  specifies  the  transformation  from  the  body  coordinate  sys¬ 
tem  to  the  earth-cer.tered  inertial  system.  The  derivation  of  the 
equations  can  be  found  in  reference  2. 


WXS  -  IVYS 
XDIF  -  YDIF 

XC11ANG  -  YCi!ANG 

XFIRST  -  XLAST 

XUW  -  YUW 
XWANT  -  YWANT 
XYZ(l)  -  XYZ(2) 


I 


Let  V  represent  the  velocity  of  the  rocket  and  let  if  represent 
the  rotation  of  the  x,  y,  z  system  with  respect  to  the  X,  Y,  2  sys¬ 
tem.  The  x,  y,  z  components  of  V  and  if  are  designated  by  u,  v,  w  and 
p,  q,  r,  respectively. 


The  equations  of  motion  derived  in  [1]  are 

mu 

a 

m(rv  -  qw)  +  F„ 

• 

mv 

s 

m(pw  -  ru)  ♦  Fy 

• 

U1W 

s 

m(qu  -  pv)  ♦  Fz 

• 

• 

!xxP 

s 

CIyy  *  qr  -  JxxP 

♦  L 

i  q 

yyH 

= 

CIzz  "  *xx>  Pr  -  V 

♦  M 

lzzT 

a 

CIxx  -  V  Pq  •  1zzT 

♦  N 

where  m  is  the  rocket  mass,  Ixx#  I  ,  are  the  moments  of  inertia, 

Fx,  Fy  and  F^  are  the  body  components  of  the  applied  forces  and  L,  M, 

N  are  the  body  components  of  the  moments.  The  forces  and  moments  are 
computed  from  the  equations 


Fx  "  Cxq'S  +T-mgx 

Fv  “  “cn  sin  0  q'  S  -  mg 

3  a  1 

Fz  "  ~Cn  sin  a*q'  s  "  m8z 

a 

L  ■  %  5  *  eg-) 

6  pa 

M  -  [Cm  sin  a*  ♦  Cm  Sd  >  m(l-cg)2q 

a  q  a 

N  «  [-C  sin  0  ♦  Cn  Sd  ♦  m(l-cg)2r 

a  q  a 


-  U 
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where 


2  J 
q’  •  1/2  Pv  is  the  dynamic  pressure, 

S  is  the  reference  area, 

d  is  the  reference  length,  and 

v  is  the  magnitude  of  the  aerodynamic  velocity. 
•  •  • 

®(l-cg)q  and  m(l-cg)r  are  jet  damping  terms. 

The  thrust,  T,  is  computed  from  the 

formula 


where 


T-T  ,  ♦  A  (P  .  -  P  ) 

s.l.  e  v  s.l.  a 


T  ,  is  sea  level  thrust 
S  •  X  • 


P  .  is  sea  level  pressure 
S»l# 


is  atmospheric  pressure 


and 


A  is  the  area  of  the  exit  nozzle, 
e 


The  body  components  of  the  acceleration  due  to  gravity  are  des 

ignated  g  ,  g  ,  and  g  .  These  are  obtained  from  the  inertial  compo 
x  y  z 

nents  of  gravity  which  are  assumed  to  be 


Y  Me  X 


[1  -  3  R2  v(l  -  3  cos  2  X)/R2] 
eq 


Y  Me  Y 

-stj- 


[1  -  3  R2  v(l  -  3  cos  2  X)/R2] 
eq 


Y  Me  Z 

-*3- 


[1  -  3  R2  (1-3  cos  2  X)/R2] 
eq 
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V 


where  y  ’! 
S 


16  3  2 

1.40775  x  10  ft  /sec 
.273  x  Iff3 


X  is  geodetic  latitude. 


Wind  enters  the  equation  in  the  terms  sin  E,  sin  a*  and  v  .  A. 

precise  formulation  is  presented  in  [1],  The  wind  components  and  the 
atmospheric  pressure  and  atmospheric  density,  p,  are  obtained  from 
a  table  using  the  subroutine  TABL.  The  wind  components  in  the  body 
system  are  obtained  by  applying  the  transformations  described  in  the 
section  "Coordinate  Systems  and  Transformations." 

The  stability  derivatives  C  ,  C  .  CB  ,  Cff  ,  C  and  C  .  which 
1  x'  n  8  lr.>  I  '  m  m  * 

a  c  p  a  q 

are  defined  in  the  section  on  Input,  are  also  obtained  by  a  table 
lookup.  The  routine  has  an  option  of  either  obtaining  C  from  the 

a 

table  or  computing  this  value  from 


The  derivatives  of  the  direction  cosines  are  computed  from  the 
formulas 


1.  = 

1 

rm. 

i 

"  n 

i  = 

1,2,3 

• 

m.  = 
i 

pn. 

-  rl. 

i 

i  = 

1,2,3 

« 

n.  = 

i 

qJi 

-  pm. 

i  = 

1,2,3 

The  inertial  component  of  the  accelerations  are  obtained  by  applying 
the  body-to-inertial  transformation. 

The  only  input  parameter  to  the  routine  which  is  not  transmit¬ 
ted  through  Common  is  the  integer  J  which  appears  in  the  calling  se¬ 
quence.  This  parameter  indicates  the  step  in  the  Runge-Kutta  integra¬ 
tion  routine  at  which  EVAL  was  called.  The  principal  flads  are: 

AE  -  Exit  area  of  motor  nozzle  (sq.  ft.). 

CC  -  Center  of  gravity  (feet  from  nose). 
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CHEK 


CLA7 

CMA 

CMQ 

CNA 

CP 

cx 

DAMP 

DMKOMK  (3) 

DSQ 

F.PSQ 

FLOW 

FM 

FX 

FY 

FZ 

G 

GCONl) 

GNU  J 

iiRO 

irr 

II! 


-  Magnitude  of  radius  vector  to  rocket  (feet). 

-  Cosine  of  geocentric  latitude. 

-  Center  of  pressure  (feet  from  nose) . 

-  Damping  moment  coefficient. 

-  Normal  force  coefficient. 

-  Center  of  pressure  (feet  from  nose) . 

-  Axial  force  coefficient. 

-  Jet  damping  term. 

-  Derivatives  of  moments  of  inertia. 

-  Reference  area  of  aerodynamics. 

-  Square  of  the  eccentricity  of  the  earth. 

-  Mass  flc,;  rate. 

-  Mass. 

-  Moment  of  inertia  about  x-axis. 

-  Moment  of  inertia  about  y-axis. 

-  Moment  of  inertia  about  z-axis. 

-  Gravity. 

-  Constants  used  in  evaluation  of  gravity  for  geodetic 
earth. 

-  Magnitude  of  radius  vector  to  racket. 

-  Height  of  rocket  above  oblate  earth. 

-  Indicates  whether  to  call  ballistic  factor  option  (2) 
or  not  (1) . 

-  F.ntry  of  wind  currently  being  used  in  ballistic  factor 
option. 


1 


I  WIND 


KAN  ••  Indicates  if  burning  (1)  or  coasting  (2)  phase. 

MAI,  -  Indicates  if  CP  (1)  or  CMA  (2)  is  input. 

NTYPF.  -  Indicates  what  ^ype  of  trajectory  being  simulated. 

p  -  x-component  of  the  rotation  of  the  body  system  with 

respect  to  the  inertial  system. 

PA  -  Atmospheric  pressure. 

PSL  -  Atmospheric  pressure  at  sea  -evel. 

q  -  y-component  of  rotation  of  body  system  with  respect 

->  inertial  system. 

QP  -  Dynamic  pressure. 

r  -  z-component  of  rotation  of  body  system  with  respect 

to  inertial  system. 

REFL  -  Reference  length  for  aerodynamics. 

REQ  -  Radius  of  earth  at  equator. 

RIIO  -  Atmospheric  density. 

RM  -  Mach  number. 

RS  -  Radius  of  oblate  spheroid  along  radius  vector  to  rochet. 

SA  -  Sine  a* 

SB  -  Sine  B. 

SLAT  -  Sine  of  geocentric  latitude  of  launcher. 

SLATG  -  Sine  of  geodetic.  latitude  of  launcher. 

T1ME0  -  Time  rocket  leaves  end  of  launcher. 

TSL  -  Thrust  at  sea  level. 

UP,  VP,  WP  -  Body  coordinate  system  components  of  aerodynamic  velocity. 
VA  -  Atmospheric  velocity. 
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VS  -  Speed  of  sound. 

w  -  Angular  rotation  of  the  earth. 

KDOT  (1-3)  -  Body  components  of  p,  q  and  r. 

•  ♦  • 

WDOT  (4-6)  -  Body  components  of  p,  q,  r. 

WX,  IVY,  WZ  -  Components  of  wind  in  inertial  coordinate  system. 

WXP  (50),  WYP(50)-  Components  of  wind  in  launcher  coordinate  system. 

XLT  -  Length  of  rocket  (feet). 


THE  INPUT  SYSTEM 


The  data  required  for  a  simulation  is  read  into  the  computer 
at  the  beginning  of  the  program.  These  data  are  stored  in  memory 
for  later  use.  The  data  required  for  a  given  phase  (a  specific  por¬ 
tion  of  the  trajectory)  are  extracted  from  this  storage  area  and  placed 
into  the  working  storage  area  at  the  beginning  of  the  phase. 

SUBROUTINE  Pintrj 

a  This  subroutine  is  used  to  read  in  all  input  data  and  to  store 
the  data  phasewise.  All  card  formats  are  in  columns  of  eight  with 
a  few  minor  exceptions.  The  first  six  columns  of  all  data  cards  can 
be  used  "or  identification  purposes.  The  first  card  of  the  input 
deck,  containing  the  number  of  distinct  trajectories  in  the  job,  is 
read  in  the  monitor  routine,  LRBM,  in  columns  17-24. 

Description  of  the  Input  Deck 

The  second  input  card  determines  the  type  of  trajectory  and  an 
indicator  to  determine  if  the  regular  printout  is  desired  or  not. 

If  column  9  contains  a: 

1  -  a  regular  trajectory  is  calculated. 

2  -  an  iteration  for  launcher  angles  is  run. 

3  -  an  angular  displacement  table  is  computed. 

4  -  the  ballistic  factors  and  unit  wind  effect  are  calculated. 

5  -  a  dispersion  analysis  is  run. 


23 


printouts 


If  also  on  the  second  card,  column  17  contains  a  1, 
occur  only  at  the  end  of  each  phase. 

The  input  for  cards  3  through  6  is  as  follows: 


Card 

Column 

Parameter 

Description  (Units) 

3 

1-72 

IDEN 

Alphanumeric  information 
used  for  identification. 

4 

9-16 

TO 

Time  (sec)  rocket  leaves 
end  of  launcher. 

17-24 

ZL 

Mean  sea  level  altitude  at 
end  of  launcher  (ft). 

25-32 

ZIM 

Mean  sea  level  altitude  of 
impact  area  (ft). 

33-40 

XLAT 

Geocentric  latitude  of  launcher 
(degrees) . 

41-48 

XLONG 

Geocentric  longitude  of  launcher 
(radians) . 

49-50 

NFAZE 

Total  number  of  phases  (right 
justified) . 

57-58 

NBST 

Phase,  end  of  which,  to  pick 
up  booster  (right  justified) . 

5 

9-16 

UO 

x  component  of  velocity  in 
body  system  at  end  of  launcher 
(ft/sec) . 

17-24 

VO 

y  component  of  velocity  in 
body  system  at  end  of  launcher 
(ft/sec) . 

25-32 

WO 

z  component  of  velocity  in 
body  system  at  end  of  launcher 
(ft/sec) . 

33 

MAL 

Indicator;  if  equals  1,  in¬ 
put  CP;  if  equals  2,  input 

CMA. 
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Card 

Column 

Parameter 

Description  (Units) 

j 

6 

9-16 

Til 

Tilt  angle  (degrees). 

17-24 

AL 

Azimuth  angle  (degrees). 

2S-32 

XWANT 

East-west  component  of  de¬ 
sired  impact  (ft). 

33-40 

YKANT 

North-south  component  of 
desired  impact  (ft). 

41-48 

DELTA 

Iteration  tolerance. 

49-54 

TOUT  I L 

T^„or  tilt  effect  (ft/radian) 

The  remainder  of  the  input  consists  mainly  of  tables.  All  tables 
use  one  card  for  each  value  of  the  independent  parameter  (height, 

Mach  number,  or  time).  The  first  card  of  each  table  contains  a  right 
justified  integer  in  columns  7  and  8  which  indicates  the  number  of 
cards  in  the  table. 

There  are  five  types  of  tables: 

1.  Atmospheric  Table 


Column 

Parameter 

Description  (Units) 

9-16 

HH 

Height  above  sea  level  (ft). 

17-24 

RHO 

Density  (slugs/ft3). 

25-32 

VS 

Speed  of  sound  (ft/sec). 

33-4C 

PA 

Pressure  (lbs/sq  ft). 

2. 

Wind  Table 

9-16 

HP 

Height  above  sea  level  (ft). 

17-24 

WXP 

East-west  component  of  wind  (i 

25-32 

WYP 

North-south  component  of  wind 
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Mach  Table 


Column 

Parameter 

Description 

j 

(Units) 

9-16 

PM 

Mach  number, 

17-24 

CX 

Axial  force 

coefficient. 

2S-32 

CNA 

Normal  force  coefficient  (per  radian). 

33-40 

CMA 

Restoring  moment  coefficient  or  center 
of  pressure  (ft  from  nose). 

41-48 

CMQ 

Pitch  damping  moment  coefficient. 

49-56 

CLD 

Roll  driving  coefficient. 

56-64 

CLP 

Roll  damping  coefficient. 

4 

.  Time  Table 

9-16 

T 

Time  (sec) 

17-24 

FX 

2 

Roll  moment  of  inertia  (slug/ft  ). 

25-32 

FY 

2 

Pitch  moment  of  inertia  (slug/ft  ). 

33-40 

TSL 

Thrust  (lbs). 

41-48 

CG 

Center  of  gravity  (ft  from  nose). 

49-56 

FM 

Mass  (slugs). 

5, 

.  Ballistic  Factor  Table 

9-16 

BF 

Ballistic  factor. 

For 

each  phase 

there  are  two  Phase  Cards  containing  in: 

1-8 

TBO 

Time  phase  terminates  (sec). 

9-16 

AE 

Area  of  motor  exit  nozzle  (sq  ft  ). 

17-24 

DSQ 

Reference  area  used  in  aerodynamic  ealeu 
lations  (sq  ft). 

25-32 

REFL 

Reference  length  used  in  aerodynamic  cal 
culations  (ft). 
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Column 

Parameter 

Description  (Units) 

33-40 

XLENTH 

Length  of  rocket  (ft)* 

| 

41-48 

eptiny) 

Integration  tolerances  for  Runge-Kutta- 

* 

f 

Gill  numerical  integration. 

i 

49-56 

EPBIG J 

I 

57-64 

FINC 

65-72 

TMII 

Maximum  integration  interval  (sec). 

| 

jil 

73-80 

PI 

Print  interval  (sec). 

1 

1-8 

TFP 

Time  of  first  print  (sec). 

£ 

9 

KAN 

Indicator;  if  equals  1,  input  Mach  and 
time  tables;  if  equals  2,  Mach  table 
only. 

V 

js: 

Depending  on  the  type  of  trajectory  desired,  several  other  cards 
may  be  necessary.  These  cards  include: 


Unit  Wind  Effects  Card 

Column 

Parameter 

Description  (Units) 

9-16 

DELC 

Cross  unit  wind  effect 

(m/mph) . 

17-24 

DELT 

Tail  unit  wind  effect 

(m/mph) . 

25-32 

DELH 

Head  unit  wind  effect  (m/mph). 

Angular  Displacement  Kind  Card 

9-16 

XFIRST 

The  first  E-W  wind  used  (mph). 

17-24 

XLAST 

The  last  E-W  wind  used  (mph). 

25-32 

XCHANG 

Increment  for  E-W  wind  change  (mph). 

33-40 

YFIRST 

The  first  N-S  wind  used  (mph). 

41-48 

YLAST 

The  last  N-S  wind  used  (mph) . 

49-56 

YCHANG 

Increment  for  N-S  wind  change  (mph) . 
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The  input  for  the  ballistic  factor  option  is  the  same  as  for  the 
regular  trajectory.  The  first  card  of  the  wind  table  should  contain 
a  zero  wind.  The  remainder  of  the  cards  of  the  wind  table  should 
contain  winds  of  constant  magnitude  and  direction. 

The  input  for  the  dispersion  analysis  is  the  same  as  for  the 
regular  trajectory  except  after  each  Mach  and  time  table,  another 
card  follows.  This  card  contains  a  percentage  of  change  in  decimal 
form  plus  one  in  the  same  field  as  the  parameter  to  be  changed.  That 
is,  if  1.1  appears  in  columns  17-24  on  the  card  following  the  Mach 
table,  the  axial  force  coefficient  is  increased  ten  percent.  If  the 
parameter  is  not  to  be  changed,  a  1.  must  appear  in  the  same  field 
as  the  parameter  itself. 

If  one  also  desires  to  bring  the  booster  to  impact,  another  phase 
must  be  added.  This  phase  follows  the  last  regular  phase  and  consists 
of : 


Phase  Cards 
Mach  Table 
Time  Table  (if  KAN 
On  card  2,  the  value  of 
SUBROUTINE  Pliasin  (J) 


Repeated  as  required. 


is  1) 


NFA2E,  should  include  the  booster  phases. 


This  subroutine  selects  the  data  required  for  the  Jth  phase  and 
places  these  data  in  working  storage  areas.  Recall  that  all  data 
input  is  at  the  beginning  of  the  trajectory.  The  principal  flads 
for  the  input  system: 


AE 

AER (10) 

CG (40) 
CGB(40, 10) 
CMA(20) 
CMAB(20, 10) 
CMQ(20) 
CM0B(20, 10) 


-  Area  of  exit  nozzle  (sq  ft). 

-  Storage  of  AE  for  all  phases. 

-  Entries  of  center  of  gravity  in  time  tables. 

-  Storage  of  CG  for  all  phases. 

-  Entries  of  center  of  pressure  in  Mach  table. 

-  Storage  of  CMA  for  all  phases. 

-  Entries  of  coefficient  of  damping  moment  in  Mach  table. 
-Storage  of  CMQ  for  all  phases. 


29 


CNA(20) 

CNA8(20,10) 


Entries  of  coefficient  of  normal  force  in  Mach  table. 


CX  C  2  0) 
CXB(20, 10) 
DSQ 

DSQB(IO) 

EPBIG  \ 
tPTINYJ 

EPBIGB (10) ) 
EPTINB(10)J 

FM(40) 

FMB(40, 10) 

FX(40) 

FXB(40, 10) 
FY (40j 

FYB(4C, 10) 
IOPTUB(IO) 
KAN 

KAN'Bf  10) 

NM 

NMB(10.t 

NT 

NTB(IO) 

PI 


-  Storage  of  CNA  for  all  phases. 

-  Entries  of  coefficient  of  axial  force  in  Mach  table. 

-  Storage  of  CX  for  all  phases. 

-  Reference  area  used  in  aerodynamic  calculations. 

-  Storage  of  DSQ  for  all  phases. 

-  Integration  tolerances  for  Runge-Kutta-Gill  numerical 
integration. 

-  Storage  of  EPBIG  and  EPTINY  for  all  phases. 

-  Entries  of  mass  in  time  tables. 

-  Storage  of  FM  for  all  phases. 

-  Entries  of  pitch  moment  of  inertia  about  x-axis  in 
time  table. 

-  Storage  of  EX  for  ail  phases. 

-  Entries  of  pitch  moment  of  inertia  about  y-axis  ;n 
time  tabic. 

-  Storage  of  FY  for  all  phases. 

-  Storage  of  IOPTUN  for  all  phases. 

-  Indicator,  if  equal  to  two,  only  Mach  table  is  nec¬ 
essary. 

-  Storage  of  KAN  for  all  phases. 

-  Number  of  values  in  Mach  table. 

-  Storage  of  NM  for  all  phases. 

-  Number  of  values  in  time  table. 

-  Storage  of  NT  for  all  phases. 

-  Prira  interval. 
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PIB(IO) 

PM (20) 
PMB(20, 10) 
T(40) 
TB(40, 10) 
TBO 

TBOB(IO) 

TFP 

TFPB(IO) 

TMI I 

■nillB(lO) 
TSL(40) 
TSLB(40, 10) 
XLT(IO) 
XLNTH 


-  Storage  of  PI  for  all  phases. 

-  Entries  of  Mach  number  in  Mach  table. 

-  Storage  of  PM  for  all  phases. 

-  Entries  of  time  in  time  table. 

-  Storage  of  T  for  all  phases. 

-  Time  at  end  of  phase. 

-  Storage  of  TBO  for  all  phases. 

-  Time  of  first  print  to  be  used  in  phase. 

-  Storage  of  TFP  for  all  phases. 

-  Maximum  integration  interval. 

-  Storage  of  maximum  integration  interval. 

-  Entries  of  thrust  in  time  tables. 

-  Storage  of  TSL  for  all  phases. 

-  Entries  of  length  of  rocket. 

-  Length  of  the  rocket. 

THE  OUTPUT  SYSTEM 


s 


f 

i 

[ 


t 
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The  output  from  the  simulation  is  a  printed  record  of  the  sim¬ 
ulation  at  discrete  preselected  times.  The  output  is  referenced  to 
the  launcher  coordinate  system.  Throe  subroutines  are  used  in  the 
output  system.  These  are:  (1)  CHKOUT  is  used  to  determine  when  out¬ 
put  is  required.  (2)  TI2L  transforms  the  data  from  the  Inertial 
System  to  the  Launcher  System.  (3)  XOUT  prepares  the  magnetic  tape 
for  the  printe?.  XOUT  and  CHKOUT  are  called  by  LRBM;  TI2L  is  called 
by  XOUT. 

SUBROUTINE  CHKOUT 

Output  from  the  trajectory  is  required  at  the  following  times: 

(1)  specified  intervals  within  a  phase;  (2)  end  of  each  phase;  and, 

(3)  impact.  The  subroutine  CHKOUT  is  called  at  the  end  of  each  in¬ 
tegration  interval  to  see  if  any  of  these  conditions  have  been  met. 


[  ' 
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An  indicator  IOUT  is  set  equal  to  2,  3,  or  4,  respectively,  if  con¬ 
dition  (1),  (2)  or  (3)  is  satisfied  and  set  equal  to  1  otherwise. 

The  logic  flow  of  LRBM  is,  of  course,  affected  by  the  value  of  IOUT. 
In  performing  the  checks  it  is  assumed  that  impact  will  not  occur 
before  a  simulated  time  of  20  seconds. 

The  integration  interval  is  also  controlled  by  C11K0UT.  This 
control  is  limited  to  control  necessary  to  preclude  "jumping  over" 

(1)  or  (2) , 

SUBROUTINE  TI2L  (KWIND,  JF). 

This  subroutine  is  used  to  convert  values  of  position,  velocity, 
and  acceleration  from  the  earth-centered  inertial  coordinate  system 
to  the  launcher  coordinate  system  for  use  in  the  subroutine,  XOUT. 

The  total  aerodynamic  velocity  and  the  azimuth  and  elevation  angles 
of  the  velocity  vector  are  also  calculated. 

The  position,  velocity,  and  acceleration  in  the  launcher  system 
are  obtained  from  the  following  equations: 

*L  "  T<*I  -  V 

*  *  -  * 

Rl  =  T(Rj  -  W  x  R  ) 

t  ii  t 

Rl  =  T(R  -  if  x  (W  x  ft  )  -  2W  x  R  ) 

where  is  the  position  vector  in  the  launcher  system, 

is  the  position  vector  in  the  inertial  system, 

W  is  the  earth's  rotation  vector, 

T  is  the  linear  transformation  from  the  inr~«.ial  to  the 
launcher  coordinate  system,  and 

R  is  the  radius  vector  from  the  inertial  to  the  launcher 
o 

system. 

KWIND  is  an  indicator  used  in  the  calling  sequence  tc  determine 
if  the  accelerations  should  be  computed  in  the  launcher  system  or 
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RETURN 


not.  Only  when  KWIND  equals  one  are  these  accelerations  computed. 

The  principal  flads  of  the  output  system  -are: 

All  -  Azimuth  angle  of  velocity  vector  (deg). 

B(l)-B(2)-B(3)  -  Launcher  coordinates  of  aerodynamic  acceleration 

(ft/sec/sec) . 

CHEK  -  Magnitude  of  radius  vector  to  rocket  (ft) . 

CLAT  -  CLATG  -  Cosine  of  geocentric  and  geodetic  latitude. 

D(l)-D(2)-D(3)  -  Launcher  coordinates  of  aerodynamic  velocity  (ft/sec). 

DIF  -  Difference  between  the  time  of  the  simulated  tra¬ 

jectory  and  the  time  of  the  beginning  of  the  phase 
(sec) . 

EMPT  -  Sum  of  the  radius  of  the  oblate  spheroid  along  ra¬ 

dius  vector  to  the  rocket  and  the  height  of  the 
impact  area  above  sea  level  (ft) . 

END PI  I  -  Difference  between  current  time  of  simulated  tra¬ 

jectory  and  time  of  the  end  of  the  phase  (sec). 

EPSQ  -  Square  of  the  eccentricity  of  the  earth. 

H  -  Integration  interval  (sec) . 

IDEN  -  Alphanumeric  information  used  for  identification. 

IJK  -  Current  line  of  print  out  on  page. 

NUF  -  Set  equal  to  1  if  CHEK  is  less  than  or  equal  to 

EMPT. 

PI  -  Print  interval  (sec). 

PO  -  If  PO  equals  zero,  print  out  occurs. 

RO  -  Radius  of  earth  (ft). 

RS  -  Radius  vector  of  oblate  spheroid  along  radius  vec¬ 

tor  to  rocket  (ft), 
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SLAT 

SLATG 

TBO 

TFp 

Til 

TIMEO 

TO 

V 

W 

WT 

WXS,  KYS 
XL,  YL,  ZL 

XYZ  (3) 

YD  (6) 

ZIM 


-  Sine  of  geocentric  latitude  of  launcher. 

-  Sine  of  geodetic  latitude  of  launcher. 

-  Time  of  end  of  current  phase  (sec) . 

-  Time  of  first  print  out  for  current  phase  (sec)., 

-  Elevation  angle  of  velocity  vector  (deg). 

-  Time  (sec)  rocket  leaves  end  of  launcher. 

-  Time  of  beginning  of  phase  (sec). 

-  Aerodynamic  velocity  in  launcher  coordinate  system. 

-  Omega,  rate  of  angular  rotation  of  earth. 

-  Distance  earth  rotates  in  time,  T. 

-  E-W  and  N-S  components  of  wind  (mph) . 

-  Translation  vectors  from  earth-centered  inertial  co¬ 
ordinate  system  to  launcher  coordinate  system. 

-  Launcher  coordinates  of  position  of  rocket  (ft) . 

-  Relative  motion  velocities  of  inertial  system  with 
respect  to  launcher  system. 

-  Impact  height  above  sea  level  of  impact  area  (ft). 

THE  INTEGRATION  ROUTINES 


The  trajectory  is  determined  by  integrating  15  differential  equa 
tions  of  motion.  Three  are  of  second  order,  while  the  other  12  are 
first  order.  Iterated  integration  is  used  for  the  second  order  equa¬ 
tions.  The  integration  is  monitored  by  LRBM  through  the  subroutine 
INTEG. 


The  subroutine  INTEG  monitors  the  mechanics  of  the  numerical 
integration  through  the  routines  RKG  and  RKERR2. 

The  SUBROUTINE  RKG  (II,  TIME)  is  used  to  integrate  the  equations 
of  motion  and  uses  the  fourth  order  Runge-Kutta-Gill  numerical  inte¬ 
gration  scheme.  RKG  calls  upon  EVAL  for  evaluation  of  the  equations 
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SAVE  INITIAL 
CONDITIONS 


RXG (H .TIME) 


NOGOOD 


RK£RR2(H) 


S  NOGOOD 


IS  H  <  .0001? 


HALF  THE 

INTEGRATION 

INTERVAL 


_|V(I,S)  *  YHALF(I) 


of  motion.  II  is  the  integration  interval  currently  being  used.  TIME 
is  the  current  tine  in  the  trajectory  simulation. 

The  SUBROUTINE  RKF.RR(ii)  is  used  to  check  the  validity  of  the 
integration.  This  check  is  performed  by  comparing  the  value,  say  Yj, 

obtained  for  an  interval  of  length  1!  with  the  value,  say  yjj»  obtained 

by  performing  two  integrations  using  the  value  11/2.  Specifically  the 
following  ratio 

(Yjj  -  is  compared  with  two  tolerance  values  and  e^. 

(This  comparison  is  made  separately  for  each  of  the  integrated  values.) 

If  (Yjj.  -  Yj)/Yjj  i-s  greater  than  for  any  of  the  parameters, 

li  is  replaced  by  U/2  and  the  integration  is  redone.  If  all  values  of 
(Yu  -  Yj)/Yjj  <  e2  t^cn  not  onlY  the  integration  accepted  but  the 

integration  interval  for  the  next  step  is  doubled. 

If  at  any  time  during  the  simulation  II  becomes  less  than  .0001, 
a  message  is  sent  to  LRBM  using  the  indicator  DONE  to  terminate  the 
simulation.  The  indicator  NOf.OOD  is  used  to  tell  whether  the  integra¬ 
tion  is  acceptable  or  not. 

SUBROUTINE  Mtrxl  (Al,  A2,  NC2,  PRO) 

This  subroutine  is  used  to  multiply  two  matrices.  Al,  A2,  NC2, 
and  PRO  are  variables  used  in  the  calling  statement.  Al  is  a  3  x  3 
matrix.  A2  is  either  a  3x1,  3x2,  or  3x3  matrix  depending  on  whether 
NC2  is  equal  to  one,  two,  or  three.  The  product  of  these  two  matrices, 
Al  and  A2,  is  then  stored  in  the  matrix,  PRO. 

SUBROUTINE  TABL  (LOOKUP,  J) 

This  subroutine  is  used  to  interpolate  linearly  for  values  from 
the  atmosphere,  wind,  Mach,  and  time  tables.  The  parameter  LOOKUP 
is  used  to  designate  which  table  is  to  be  used. 

To  conserve  computer  time  the  routine  is  programmed  to  "remember" 
the  arguments  used  in  the  previous  entry  to  the  table.  The  "search" 
is  either  forward  or  backward  depending  upon  the  present  value  of  the 
independent  variable. 
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Suppose  the  independent  variable  lies  outside  the  arguments  of  the' 
table.  For  values  from  the  Mach  or  atmosphere  tables  linear  extra¬ 
polation  is  used;  for  the  wind  tables  the  value  is  set  to  zero;  and, 
for  the  time  tables  the  last  value  of  the  table  is  chosen. 

The  interpolation  routine  for  the  time  table  also  provides  for 
the  calculations  of  the  derivatives  of  the  moments  of  inertia,  and 
of  the  mass  flow  rate.  The  principal  flads  are: 

CD  -  Interpolated  value  of  axial  force  coefficient. 

CG  (40)  -  Entries  of  center  of  gravity  in  time  table. 

CL  -  Interpolated  value  of  normal  force  coefficient. 

CLD  -  Roll  driving  moment. 

CLP  -  Roll  damping  moment. 

CMA(20)  -  Entries  of  center  of  pressure  in  Mach  table. 

CMQ(20)  -  Entries  of  coefficient  of  pitch  damping  moment  in 

Mach  table. 

CMT  -  Interpolated  value  of  coefficient  of  pitch  damping  mo¬ 

ment. 

CNA(20)  -  Entries  of  normal  force  coefficient  in  Mach  table. 

CND(40)  -  Entries  of  roll  driving  coefficient. 

CNP(40)  -  Entries  of  roll  damping  coefficient. 

CON  -  Constant  used  for  interpolation. 

CP  -  Interpolated  value  of  center  of  pressure. 

CX(20)  -  Entries  of  coefficient  of  drag  in  Mach  table. 

DENS  -  Interpolated  value  of  atmospheric  density. 

_pDMK0MK(3)  -  Derivatives  of  moments  of  inertia. 

FM(40)  -  Entries  of  mass  in  time  table. 

FX (40)  -  Entries  of  moment  of  inertia  about  x-axis  in  time 

table. 
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FY (40) 
ltt!(44) 
IIP(SO) 

HT 

11,12,13,14 

LOOKUP 

N 

NM 

NT 

NW 

PA (44) 

PM  (20) 

PRES 

R1K>(44) 

RM 

SPD 

T(40) 

TME(6) 


TSL(40) 

WXP(50) 

WXS 

WYP(SO) 

WYS 


Entries  of  moment  of  inertia  abput  y-axis  in  time  table. 
Entries  of  height  in  atmospheric  table. 

Entries  of  height  in  wind  table. 

Height  of  rocket  above  oblate  earth. 

Entry  last  used  in  the  various  tables. 

Indicates  which  table  to  interpolate. 

Number  of  entries  in  the  atmospheric  table. 

Number  of  entries  in  the  Mach  table. 

Number  of  entries  in  the  time  table. 

Number  of  entries  in  the  wind  table. 

Entries  of  pressure  in  atmospheric  table. 

Entries  of  Mach  number  in  Mach  table. 

Interpolated  value  of  pressure. 

Entries  of  atmospheric  density  in  atmospheric  table. 

Macli  number. 

Interpolated  value  of  speed  of  sound. 

Entries  of  time  in  time  table. 

Interpolated  values  of  moment  of  inertia  about  the  x- 
and  y-axes,  thrust,  center  of  gravity,  mass,  and  center 
of  gravity  of  the  propellant. 

Entries  of  thrust  in  time  table. 

Entries  of  E-W  components  of  wind. 

Interpolated  value  of  E-W  components  of  wind. 

Entries  of  N-S  components  of  wind. 

Interpolated  value  of  N-S  components  of  wind. 


( 


i 


40 


141 


LISTING  OF  PROGRAM  INSTRUCTIONS 


SATuNO  JOGOO. 77777.fe.ouMP 

^EXECUTE  ID JOB 

4  I 3 JOB  GO  « LOO  I  C  .  MAP 

MRFTC  LftHM  LI  ST  .PEP. DECK 

COMMON  Y (22*b  > .DY  <22.5 ) , ROFF1  ( 22 . ^ ) 

C 

C  IN  ALL  COMMON  STATEMENTS  l)  =  OW  I  VtH  .  A  =  AOPT  UN  •  8  =  lJFACT  ♦  E=cV  AL  •  C  =  CHK  *  1*  INT 
C  R2=RKEWR2»H=RKG.T=TAcJL.K*PHASlN.PsPlNTRj.X*X0ur .AND  l=T I 2LP 
C 

COMMON/OA/  IXwO.  I  Yji/D.SI  THP.COTMP.ENF  «  THP.  ALP 

c 

CCMMON/DAP/TH.AL.XdANT  * y*ant.delc.delt .delh. towt il 

r 

common/oaep/ntype 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


COMMON/OAETX/  tfXS.tf YS 

COMMON/DB/  YST0W(22  )  .OYSTOR  (22  >  .RJFSTW<22  )  «  I  tiFENO 

common/dbf/i wind* ih.bhi 

COMMON/OB TP/NW . HP ( 50  > 

C0MM0N/0£/REQ*PSL .gconi .Gnu 

CO.MMON/DEL/  SLAT.CLAT  .^lATG.ClATG. T  lMbOt« 

COMMON/OCP/  TO 

COMMON/DC/  NUF.IOUT 

COMMON /DC F/P I .TOO 

COMMON/OCI/  H 

COMMON/O I /  DONE 

COMMON/D  IX/  TIME 

COMMON/D  I R2/  NOGOOD 

COMMON/DT/  11.12.13.14 

COMMON/OP/  UO.VO. *0.2L.XLAT .NFAZE.NBST. IFOuT 

COMMON/DF/  TMIl 

COMMON/DX/  1JK 

COMMON/DL/  PO 

COMMON/DEC /EPSO 


00020 
00030 
00040 
OOObO 
00060 
00070 
00060 
00090 
00  100 
001  10 
00120 
00130 
00140 
00  ISO 
00160 
00170 

ooiao 

00190 

00200 

00210 

00220 

00230 

00240 

00250 

00260 

00270 

00280 

00290 

00300 

00310 

00320 

00330 

00340 

00350 

00360 

00370 

00380 

00390 

00400 

00410 

00420 

00430 

00440 

00450 

00460 

00470 

00400 

00490 

00500 

00510 

00520 

00530 

00540 
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DIMENSION  AICI3.3) .B1CI3.3) 

C 

C  INITIALIZE 

C 

CNTR=0. 

READ ( 5 « 600  1  XNO 
600  FORMAT! I 6XF8.0 ) 

101  CALL  PINTRJ 

DATA  REQ.PSL. tf « EPSQ . GCON 1 . GNU/2092642a . .2116. . .72921 E-04. .00672207 
1.1.40775E  16..273E-03/ 

TH=TH*.Ol 745329 
AL=AL*. 01 745329 
XLAT=XLAT  *.01745329 
SLAT=SIN<XLAT1 
CL AT «C0S  <  XLAT I 
EPSQ=EPS0**2 
TEN=EPSQ*SLAT  *CL AT 
F I VE= I . -EPSQ* I  CL AT**2 I 
DLTA  =  ATAN2 (TEN. FIVE  I 
XLATG«XLAT+DLTA 
SLATG=S1N(XLATG1 
CL ATG=COS ( XLATG ] 

ROARED/ ( 1 ,+EPSO*SLAT*SLAT/< 1 .-EPSQ) I**. 5 
ROA=RO+ZL 

Bicn  .  i  i*-i  .o 

B I C ( 2  « 1  1=0.0 
8103.1 1=0.0 
B I C  < 1  .21=0.0 
B I C (2 .2  I =— SLATG 
BIC(3«21=CLATG 
BIC  < 1 .3  1=0.0 
BIC!2.3l=CLATG 
BIC !3. 31= SLATG 
T I MEO=TO 

SET  UP  FIRST  APPROXIMATION  FOR  ITERATION  FOR  LAUNCHER  SETTINGS 

IF  (NTYPE.EQ.21  GO  TO  105 
IF  (NTYPE.NE.3I  GO  TO  1 
I YWD»1 
!XWD=I 
THP»TH 
ALP* AL 

1  05  RANGE  =  SORT  (  XWANT *XJ*ANT  +YWANT* YWANT  1 
S I THP«XWANT2RANGE 
COTHP-YWANT/RANGE 
CALL  AOPTUN  < I  1 
WR=9Y3*C0THP+WXS*S I THP 
IF  (WR.GE.O. 1  GO  TO  102 
OETH=DELT 
GO  TO  103 

102  DETH*DELH 

103  YUW=DETH*COTHP-DELC*SITHP 


00550 

00560 

00570 

00580 

00590 

00600 

00610 

00620 

00630 

00640 

00650 

00660 

00670 

00680 

00690 

00700 

00710 

00720 

00730 

00740 

00750 

00760 

00770 

00780 

00790 

00800 

ooaio 

00820 

00830 

00840 

00850 

00860 

00870 

00880 

00890 

00900 

00910 

00920 

00930 

00940 

00950 

00960 

00970 

00980 

00990 

01000 

01010 

01020 

01030 

01040 

01050 

01060 

01070 
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uu  u  o  o 


xuw*oetm*sithp+dq.c*cothp  oioao 

X0IF«*XS*XUW*528  .  01090 

YD1F»¥YS*YUW*328  •  01100 

STCA»SINITH  1*C0S  1  AL  )  OHIO 

COST*COS(TH)  01120 

TH1«ATAN2CSTCA.C0ST)  01130 

STSA«S1NCTH>*SIN(AL>  01140 

TH2*ATAN2  I  STS A . COST  1  01150 

TH1 «TH1-Y0IF/T0WTIL  01160 

TH2«TH2-XDIF/T0WTIL  011?0 

STH2*SIN  I TH2 1  01180 

CTH1«C0SITHI )  01190 

CTH2*C0S(TH2)  01200 

STH1-SINITH1 »  01210 

STCT21 «STH2*CTM1  01220 

CTS21 »CTH2*STH1  01230 

AL*ATAN2 ( STCT2 1 « CTS21 )  01240 

COAST  «COS ( AL )  01250 

STOCT 1 *STH1 /CTH1  01260 

BSTOR»A8S  ( STOCT 1 /COAST >  01270 

TH«ATAN<BSTORl  01280 

01290 

CNTR. WHICH  TRAJECTORY  IN NUF= 1  .  IMPACT 1 JK  IS  LINE  OF  OUTPUT  ON  PA  01300 

01  310 

CONTINUE  INITIALIZING  01320 

01330 

100  NUF-0  01340 

IJK»0  01350 

I  1 « 1  01 360 

I2«l  01370 

I 3« I  01380 

I  4  » 1  01  390 

Y { 1 ,5)=TIMEO  01400 

STM*S1NITH)  01440 

CTH-COSITH)  01450 

SAL-SINIAL)  01460 

CAL*C0S(AL)  01470 

Y«17.5)»0.0  01480 

Y« 18.5)*ROA*CLAT  01490 

Y< 19.5>=R0A*SLAT  01500 

AICI 1 « l 1=STH*SAL  01510 

AIC (2. 1 »»STH*CAL  01520 

A  ICO*  I  >=CTH  01530 

A  1 C 1 1  *  2 1 =  CAL  01540 

AIC (2.2)*- SAL  01550 

A l C I  3  *  2 1 *0*0  01560 

AICI 1 .3)»SAL*CTH  01570 

AIC(2.3)»CAL*CTH  01580 

AICI3.3I—STH  01590 

CALL  MTRXL  <BIC.MlC«3.Yia.5>)  01600 

OYI I  7.5 )=Y 18.5) *UO-W*Y( 18.5 ) 

0Y<18.5)=Y(9,5>*U0+W*Y(l7.5) 

DYI 1 9.5) *Y( 10.5) *U0 
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noon  non  non  non  on  nnn 


/ 


Y(2.5)=V(8,5)*OYU7.5)+r(9.5)*DYU8.5HYI  1 0  .  5  1  *DY  <  1  9 .5  ) 
Y13.5)=Y<1 1 .5)*0Y<17.5)+Y< 12.5 )*DY< 18.5) +YC 13.5 )«DY< 19.5) 
Y (4.51 =  Y 11 4.5 )*0Y< 17.5 )+Y( 15.5 >*0Y< IB. 5 >  +  Yl 16.8 l*OY< 19.5) 
Y<5.5>=0. 

Y<6.5)*0. 

Y 1 7 .5 ) =  W 
Y  <20.51 *0Y (17.5) 

Y ( 2 1 .5)=0Y(18.5> 

Y ( 22  *5  >  =  DY 1 19.5) 


DO  3  1=2.22  01640 

ROFFI ( I ,51=0.0  01650 

3  CONTINUE  01660 

CALL  PHAStN(l)  01670 

OO  2000  1=1.22  01680 

YSTOR (I  1  =  Y 1  I , 5 )  01690 

D YSTOR ( I ) *DY ( 1.5)  01700 

ROFSTR  < I  1 =R0FF 1(1.51  01710 

2000  CONTINUE  01720 

Jl»l  01730 

IF  (NTYPE.NE.4)  GO  TO  1  01740 

01750 

I W 1  NO  WHAT  ENTRY  AT  IN  BALL  FAC  SUBROUTINE  (HEIGHT  TABLE)  01760 

01770 

IWINO-1  01780 

32=1  01790 

J1  IS  THE  PHASE  WE  ARE  CURRENTLY  WORKING  IN  01800 

J2  IS  THE  PHASE  (BF)  TO  PICK  UP  TRAJECTORY  AT  FOR  8F  PROGRAM  01810 

11  BHI=HP( IWINOI  01820 

01830 

ih«2  in  ballistic  factor  routine  oie4o 

01850 

10  1H«1  01860 

CALL  EVAL151  01870 

CALL  XOUTC5.0)  01880 

01890 

MONITOR  THE  INTEGRATION  01900 

01910 

1  H.AMINllPI.  TM 11)  01920 

SO  IF  (1H.NE.2)  GO  TO  990  01930 

2  CALL  BALFACU)  01940 

IFdH.EQ.il  TMI  |»TM!  1*10.  0194E 

J2-J1  01950 

990  T 1  ME  *  Y 1 1 , 5 1  01960 

01970 

N0G000=0  INTEGRATION  INTERVAL  ACCEPTED  >1  SET  H=H/2.  01980 

01990 

99 1  NOGOOO=0  02000 

5  CALL  INTEG  02010 

IF  ( DONE . EQ . 2 . 1  GO  TO  9  02020 

02030 

MONITOR  OUTPUT  02040 

02050 

I  OUT* 1 »  NO  OUTPUT — I  OUT =2.  PRINT  OUTPUT  02060 


\ 


I 

i 
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10uT=3.  ENO  OF  PHASE — 10UT*4.  IMPACT  02070 

02080 

CALL  CHKOUT  02090 

GO  TO  (50.6.7.8). I0UT  02100 

9  CALL  XOUT  IS.  I  I  02110 

GO  TO  B  02120 

02130 

IFOUT-1 .PRINT  OUT  ONLY  AT  END  OF  PHASE  02140 

02150 

6  IF  IIFOUT.EO.il  GO  TO  5  02160 

CALL  XOUT (5.0)  02170 

GO  TO  50  02180 

7  JlaJl+1  02190 

CALL  XOUT (5.0  I  02200 

02210 

NFAZE  IS  THE  NUMBER  OF  PHASES  USED  (TOTAL!  02220 

N0ST  PHASE  END  OF  WHICH  BOOSTER  IS  TO  BE  PICKED  AT  02230 

02240 

IF  I J1 .GT. NFAZE I  GO  TO  1100  02250 

IF  (N0ST.EQ. ( Jl-1 » >GO  TO  20  02260 

25  CALL  PHAS IN I 32 1  02270 

CALL  EVAL(5>  02280 

CALL  XOUT  <5.01  02290 

GO  TO  1  02300 

02310 

STORE  V  AND  DV  ARRAYS  FOR  BOOSTER  TO  GROUND  02320 

02330 

20  DO  160  K= 1 , 22  02340 

YST0R(KI*Y(K.5»  02350 

160  CONTINUE  02360 

OO  161  L=2.22  02370 

DYSTORIL I *0Y(L.5>  02380 

ROFSTRIL I »ROFFl <L.5 I  02390 

161  CONTINUE  02400 

GO  TO  25  02410 

8  CALL  XOUT 1 5.0  I  02420 

IF  CNTYPE.EQ.4I  GO  TO  12  02430 

IF  INTYPE. EO. 31  GO  TO  13  02440 

IF  INTYPE. EO. 21  GO  TO  13  02450 

IF  INBST.E0.01  GO  T3  11  02460 

C  UNSTORE  Y.OY.T  ETC  FOR  BOOSTER  02470 

C  02480 

DO  151  K=1.22  02490 

151  Y(K.5)«YST0R(K)  02500 

DO  152  L»2«22  02510 

0YIL.51-DYST0RIL1  02520 

ROFF1  1L.51«R0FSTR(I..  1  02530 

152  CONTINUE  02540 

J1=J|41  02550 

TO*Y 11.5)  02560 

CALL  PHAS IN  I J 1  1  02570 

CALL  XOUT I  5.01  02580 

GO  TO  1  02590 
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13  CALL  A0PTUN<3)  02600 

02610 

ENF-1  FINISHED  02620 

02630 

IF  (ENF.EQ.l.l  GO  TO  1 1  02640 

GO  TO  100  02650 

02660 

IH  IS  INDICATOR  TO  CHECK  IF  WE  ARE  PAST  HEIGHT  OF  CURRENT  BF  *02670 

02680 

12  CALL  BALE AC (2  I  02690 

JI*J2  02700 

02710 

I  WIND  IS  THE  ENTRY  OF  WIND  TABLE  CURRENTLY  USING  02720 

02730 

IF  < IWIND.GT.NW1 IBFEND-1  02740 

IF  <IBFEN0.EQ.11  GO  TO  1010  02750 

CALL  PHAS IN ( J 1  I  02760 

TMI I*TMI 1/10.  02761 

H-AMIN1  (PI. TMI I)  02762 

TO-TIMEO  02763 

CALL  CHKOUT  02764 

11*1  02770 

12*1  02780 

13*1  02790 

14-1  02800 

!H*2  02810 

BH1 *HP< IWIND1  02820 

CALL  EVAL<51  02830 

CALL  XOUT  <  5. 0  I  02840 

GO  TO  990  02850 

1010  CALL  BALFACd!  02860 

1100  CNTR—CNTR+l .  02870 

IF (CNTR.LT .XNO I  GO  TO  101  02880 

RETURN  02890 

END  02900 

4IBFTC  AOPSUN  LIST. REF. DECK  02910 

SUBROUTINE  AOPTUNIKl  02920 

C  02930 

C  IN  ALL  COMMON  STATEMENTS  D*DR I VER . A* AOPTUN . B*BF ACT . E-EVAL « C»CHK . I  *  I NT  02940 

C  R2*RKERR2.R*RKG.T  =  TABL.F*PhASIN.P*P1NTRJ.X»X0UT. AND  L* T 1 2LP  02950 

C  02960 

COMMON/OA/  IXWD. IYWD.SlTHP.COTHP.ENF.THP, ALP  02970 

C  02980 

COMMON/DAP/TH , AL . XW ANT  «  YWANT . DELC . CELT . QELH, TOWT I L  02990 

C  03000 

COMMON/OAEP/NTYPE  03010 

C  '03020 

COMMON/O AETX/  WXS.WYS  03030 

C  03040 

COHMON/AP/  XFIRST.  F IRST • XCHANG. YCHANG . XL AST . YLAST . N8F ,8F < 50 > ,  03050 

1 DELTA  03060 

C  03070 

COMMON/ ABETP/  WXP ( 50 ) . WYP <50 1  03080 
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C  03090 

COMMON/ ABXL/-  XYZ(3>  03100 

C  03110 

DIMENSION  THETA ( I  I  ,  I  1  I .ALPHA < 1  I .  I  1  )  03120 

GO  TO  ( 1 .2.3) ,K  03130 

1  IF  (NTYPE.E0.2I  GO  TO  6  03140 

03150 

SET  UP  INITIAL  HIND  PROFILE  03160 

0317c 

HXS-XMRST  03180 

HYS«YFIRST  03190 

RETURN  03200 

03210 

SET  UP  SUBSEQUENT  HIND  PROFILE  03220 

03230 

2  IF  (HXS.E0.XLAST1  GO  TO  10  03240 

IXHD-lXHO+1  03250 

HXSaHXS+XCHANG  032o0 

GO  TO  21  03270 

10  IF  (HYS.GE.YLAST I  GO  TO  II  0328o 

IYWD”IYWD*1  03290 

IXHO-1  03300 

HXS-XFIRST  03310 

HYS-HYS+YCHANC  03320 

GO  TO  21  03330 

03340 

OBTAIN  BALLISTIC  HINDS  FOR  ITERATIVE  03350 

03360 

60  SUMX>0.  03370 

SUMY«0.  03380 

DO  61  I *1 , NBF  03390 

SUMX-SUMX+HXPI  1  )*BF( I  I  03400 

SUMY«SUMY+WYP(  I  1*8F(  I  1  03410 

61  CONTINUE  03420 

HXSaSUMX  03430 

HYS-SUMY  03440 

RETURN  03450 

03460 

PRINT  OUT  ANGLES  03470 

03480 

11  DO  75  M I  *  1 » I YHD  03490 

DO  75  Nlcl.lXHD  03500 

THETACNl .Ml  1  =  THETA(N1 .Ml  1*57.29578  03510 

ALPHA ( N l .Ml  I  a ALPHA ( N 1 .Ml  1*57.29578  03520 

75  CONTINUE  03530 

DO  71  M.l.IYWD  03540 

WRITE  (6.70)  (THETA ( l .M> . ALPHA ( 1 ,M) . I = | , IXWD1  03550 

70  FORMAT  ( 1H0.2X. I  I  (F8.3.2X )✓ ( I  1  (2X.F8.3  1  )  >  03560 

71  CONTINUE  03570 

03580 

ENF=1  FINISHED  TABLE  03590 

03600 

12  ENF« 1 .  03610 
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RETURN  03620 

03630 

CHECK  IF  IMPACT  IS  WITHIN  TOLERANCE  03640 

03650 
03660 
03670 
03680 
03690 
03700 

ITERATIVE  TECHNIQUE  FOR  COMPUTING  LAUNCHER  SETTINGS  03710 

03720 
03730 
03740 
03750 
03760 
03770 
03780 
03790 
03800 
03810 
03820 
03830 
03840 
03850 
03860 
03870 
03880 
03890 
03900 
03910 
03920 
03930 


STORE  LAUNCHER  SETTING  FOR  ANGULAR  OUTPUT  IN  TABLE  03940 

03950 

THETA ( IXWD. IYWDI-TH  03960 

ALPHA < IXWD. IYWD1-AL  03970 

GO  TO  2  03980 

03990 

DETERMINE  APPROXIMATE  THEORETICAL  DISPLACEMENT  FOR  NEXT  TRAJECTORY  04000 

C  04010 

21  IF  ( WXS.EQ.O. • AND. VYS - Eq • 0 • I  GO  TO  23  04020 

IF  (IYWO.EQ.il  GO  TO  22  04030 

TH- THETA! IXWD. I YWD-1 1  04040 

AL-ALPHAC IXWD. IYWD-1 1  04050 

IF(COTHP.GE.O. 1  GO  TO  41  04060 

DETH-OELT  04070 

GO  TO  42  04080 

23  THETAC IXWD. IYW01-THP  04090 

ALPHA ( IXWD. I YWD1»ALP  04100 

GO  TO  2  04  I  10 

41  DETH=DELH  04120 

42  YUW=DELT*COTHP-DELC*SITHP  04130 

YDIF-YCHANG*YUW*5280.  04140 


30  STCA«SIN(TH)*COS(ALl 
COST-COS (TH1 
THI «ATAN2<STCA.C0ST 1 
STSA-SINI THI*SIN(AL I 
TH2 -ATAN2 ( STS A .COST  1 
TH! -Thl-YOIF/TOWTIL 
TH2-TH2-XD I F/TOWT I L 
STH2-SINITH21 
CTHl -COS (THI I 
CTH2-C0SITH2) 

STH1— SINCTH1 1 
STCT21 »STH2*CTH1 
CTS2I =CTH2*STH1 
AL-AT  AN2 ( STCT2 1 . CTS2 1  1 
COAST-COS <AL> 

ST0CT1 -STH1XCTHI 

BSTOR- ABS  ( STOCT l XCOAST 1 

TH-ATANIB5TORI 

RETURN 

20  IF  (NTYPE.EQ.21  GO  TO  12 


1  XOIF-XYZI I 1-XWANT 
YD IF-XYZ 1 2  I -YWANT 

TEST VA -SORT  <  XO I F *X0  I F  »■  YO I F*  YD  1 F  > 
IF  (TESTVA.LE  DELTA)  GO  TO  20 


i 

1 


4 

! 
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XOIF«0.  04150 

GO  TC  30  04160 

22  TH»TH£TA<  IXWO-l  .  1YW01  04170 

ALaALPHAf IXWO-t . IYWD)  04180 

IF :S1THP.GE.0. )  GO  TO  Si  04190 

DETH=D£LT  04200 

GO  TO  52  04210 

51  DETH-OELH  04220 

52  XUW«DELC*COTHP+DETH*SITHP  04230 

Y0IF=0.  04240 

XDIF=XCHANG*XUW*5280.  04250 

GO  TO  30  04260 

ENO  04270 

*1BFTC  BFACT  LIST «REF .DECK  04280 

subroutine  balfaciintai  04290 

COMPTON  YC 22.5) .OY(22.Sl .ROFFI (22.51  04300 

C  04310 

C  IN  ALL  COMMON  STATEMENTS  D«DR I VER . A*AOPTUN.B«0FACT .E=EV AL . C«CHK . 1 » I NT  04320 
C  R2»RKERR2.R=RKG.T»rA8L.F=PHASIN.P*PINTHj.X«X0UT,AN0  L*TI2LP  04330 

C  04340 

COMMON/DBr  YSTORC221 .DYST0RI22 I .ROFSTR(22 > . 1BFEND  04350 

C  04360 

common/obe/iwino. im.bhi  04370 

C  04380 

C0MMON/DBTP/NW.HP(50>  04390 

c  04400 

COMMONSABETP/  WXPI50I .WVPfSOl  044 10 

c  04420 

COMMON/’ABXL/  XYZ«3>  04430 

C  04440 

COMMON/BET/  HT  04450 

C  04*60 

COMMON/B 1 R2X  YSAVEC221.DYSAVEI221 .R0FFSC221  044  TO 

C  04480 

COMMON/BPX/  IOEN  0*490 

C  04500 

O I  MENS  I  ON  XIMPAT(50|.YIMPAT(50)  04510 

C  04520 

C  BALFAC  MONITORS  BALLISTIC  FACTOR  PROGRAM  04530 

C  IBFEN0=1  WE  FINISHED  LAST  TRAJECTORY  04540 

C  04550 

IF  ( IBFENO.EO.l )  GO  TO  1 02  04560 

IF  IINTA.EQ.2I  GO  TO  15  04570 

C  04580 

C  CHECK  HEIGHT  OF  SIMULATED  TRAJECTORY  VERSUS  HEIGHT  OF  CURRENT  WIND  PR  04590 

C  04600 

IF  SHT.GE.BH1 |  GO  TO  18  04610 

GO  TO  I  0 1  04620 

C  04630 

C  STORE  Y  AND  DY  ARRAYS  FOR  NEXT  TRAJECTORY  0464Q 

C  04650 

180  DO  103  1=1.22  04660 

103  YSTORf I  >  =  YSAVE<  I  )  04670 


i 

i 


SO 


*Wf*r^r* :  wwr  ^WJPW; 


c 

c 

c 


c 

c 

c 

c 

c 


DO  104  M=2.22 
D YSTOR <  M ) -DVSAVE  (  M 1 
ROFSTR ( M I *R0FFS (  M I 

104  CONTINUE 
I  H*  1 

101  RETURN 

calculate  ballistic  factors  and  unit  wind  effect 

102  WRITE  (6.110)  IDEN 

difx*ximpat(nw>-ximpat(i i 

DIFY*Y1MPAT(NW)-YIMPAT(1 1 
DIFX*01FX/5280. 

DtFY*DIFY/52S0. 

DELR-SORT ( D IFX**2+0 I FY**2 ) 

WIND* SORT (WXP(NW)**2+WYP(NW )**2 > 

UN  I T  ■DELR/'W  I ND 
DENOM*0. 

00  105  I*1.NW 
D I FX*X I MPAT  < 1  )— X1MPAT  < 1  > 

DIFV«yimPAT<I I-Y1MPATU ) 

DIFX-DIFX^S280. 

D I FY*D1 FY/S280. 

RANGE*SORT(OIFX**2+DIFY**2 ) 

RAT  I O-RANGE/DELR 
BFACT  «RAT I O-OENOM 

WRITE  (6.111)  HP (I  ) . X I MPAT  < I ) , Y I MPAT (I  ) . D I FX • D! FY .RANGE .RAT  10. 
IBFACT 
OENOM-RATIO 

105  CONTINUE 

WRITE  (6.120)  UNIT 

120  FORMAT  (IHO/’IHOI  X17HUNIT  WIND  EFFECT =F 1 0 . 6 > 

RETURN 

STORE  CURRENT  IMPACT  AND  UNSTORE  Y  AND  DY  ARRAYS  FOR  NEXT  TRAJECTORY 
IWIND  IS  THE  ENTRY  OF  WIND  TABLE  CURRENTLY  USING 


150  1  WIND* I W I ND  +  1 

X  I MPAT ( 1 W I NO— 1  1*XY2<1  I 
Y I MPAT ( 1  WIND-1  ) *XYZ (2 ) 

DO  151  K*  1 .22 
Y (K  .3 ) *YSTOR( K ) 

1 5 1  CONT I NUE 

DO  152  L*2.22 
0Y(L.5)*DYST0R(L) 

R0FF1 (L.5)*R0FSTR(L) 

152  CONTINUE 
GO  TO  101 

110  FORMAT  ( l HI /l HOI 2A6/1 H0/1H0 10X . 6HHE I GHT  « 1 3X . I HX* 13X. IHY. 14X. 
lOELTA  X . BX . 7H0ELT A  Y. 9x .5HRANGE • 1  OX • 5HRAT 10. 1 1 X . 4HB .F./ I  HO > 

111  FORMAT  ( 1H05XF10.2.7F15.4) 

END 

SI 


7h 


04680 

04690 

04  700 

04710 

04720 

04730 

04  740 

J4  7S0 

04760 

04770 

04  780 

04790 

04800 

04810 

04820 

04830 

04840 

04850 

04860 

04870 

04880 

04890 

04900 

04910 

04920 

04930 

04940 

04950 

04960 

04970 

04980 

04990 

05000 

05010 

05020 

05030 

05040 

05050 

05060 

05070 

05080 

05090 

05100 

05110 

05129 

05130 

05140 

05150 

05160 

05170 

05180 

05190 

05200 
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*IBFTC  EVAlU  LI  ST. REE. DECK  05210 

SUBROUTINE  EVAL  (JJ  05220 

COMMON  V (22 .5 ) .DY<22.5( .R0FF1 (22.5 >  05230 

C  05240 

C  IN  ACL  COMMON  STATEMENTS  D« DR  I VER . A» AOPTUN . B »8F ACT . E »EV AL . C"CMK . I « I  NT  05250 
C  R2*RKERR2.R=RKG. T« TASL .F *PhAS I N .P«P 1 NTR J. X*XOuT . AND  L»1I2LP  05260 

C  05270 

COMMON/DAEP/NTVPE  05280 

C  05290 

C0MM0N/0AETX2  WXS.WYS  05300 

C  05310 

COMMON/'OBE/'IWIND.  IH.BHI  05320 

C  05330 

COMMON/DE/REO.PSL.GCONI .GNU  05340 

C  05350 

COMMON20EL/  SLAT . CL AT . SLATG . CL ATG • T I MEO. W 
C  05370 

COMMON/AHETP/  WXPC50 1 .WYPC50I  05300 

C  05390 

COMMON/BET/  hT  05400 

C  05410 

COMMON/EC/  CHEK.RS  05420 

C  05430 

COMMON/EP/  MAL  05440 

C  05450 

COMMON/OEC2EPSQ  05460 

C  05470 

COMMON/ET /WHO .VS.PA.CX. CNA . CMA . CMO. Fx  «  FY . TSL . CG • FM . RM . DMKOMK (31 . 

iflo^.cld.clp 

C  05500 

COMMON/EFT/KAN  05510 

C  05520 

COHMON/EE/  REFl. AE.DSQ.XLT.FINC 

05540 

DIMENSION  A (3.3 ) .F (3) »E(3) . AMI3.3) .WS <3  I  05550 

EQUIVALENCE  (E (  I  I . WX ) • (E (2  I >WY I . (E ( 3 ) . WZ )  05570 

0  1  MENS  I  ON  WOOTI61 

EQU1 VALENC  ( WDOT (1  I .P ) ,  IWDOT  <2 ) .0  I . (HOOT (3 ) .R I 

05580 

NORMALIZE  DIRECTION  COSINES  05590 

C  05600 

DO  1  1=8.14.3  05610 

UP =1+2  05620 

AY. SORT ( Y ( I  . J]*»2  +  Y (  I  +  1  ,  JI*»2  +  Y( JP. J>**2  >  05630 

DO  1  K  = I • JP  05640 

l  Y(K. J >=YIK. J)/AY  05650 

C  05660 

C  INTERPOLATE  FOR  4  I  NOS  AND  COMPUTE  COMPONENTS  IN  BODY  SYSTEM  05670 

C  05680 

CHEK *  SORT ( Y ( l 7. J  >*»2+Y< 18. J  >**2  +  Y( 1 9. J  >*»2 I  05690 

RS  =  20926428, / (  ( 1 .+EPSQ*<  Y< 1 9. J I/CHEK )*»22 ( 1 .-EPSCI) >*».5>  05700 

MRO  =  SORT ( Y < 1 7. J >»*2+Y ( 1 8. J l*«2+Y( 19, J>**2 >  05710 

HT=HRO-RS  05720 
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CALL  TABL < 1 « J  > 

IF  <NTVPE.E0.3»  GO  TO  12 
IF  <NTVPE.NE.*>  GO  TO  II 
IFIIH.NE.2I  GO  TO  13 
IF  (HT.GT.8MI I  GO  TO  13 
WXS*WXP< IWIND) 

WYS-W YP<  I  W 1  NO ) 

GO  TO  12 

I  ;l  WXS-O. 

W YS»0 • 

GO  TO  12 

II  CALL  TABL  <2. J ) 

12  WS< I )«FX S 

WS<2I*WYS 

WS<3>*0. 

WT«W*<Y<1 .JI-TIMEOI 

COMPUTATION  OF  LAUNCHERTO  EC  I  ROTATION  MATRIX 


A  < 1 . 1  ) ~  —COS  <  WT I 
A<2. 1  >=  — S INC  WT ) 

A  (  3 . I  1=  0.0 
A  <  3 .2  Is  CLATG 
A <3.31 »  SLATG 
A< 1 .2 )>  — A<3.3)*A<2. 1  I 
A  <2.2 1“  A  <  3 • 3  I  *A  < 1 .1  I 
A  < I • 3 ) *  A<3.2I*A<2.1  I 

AC2.3I*  —A  <3.21 *A  <1.1  I 
CALL  MTHXL  CA.WS.l.F) 

F< 1 )*-F< 1 1*1 .466667-W*Y< 18. Jl 
F  <2  )*— F 12  1*1 .466667+W*Y  < I  7, Jl 
F<3I»-F<3I*1 .466667 
WX 1  NT  *F  < 1  I 
W Y I NT  *F  <  2  I 
WZ I NT *F  <  3 ) 

DO  2  1*1.3 

AMI  1 . I )*Y< 1+7. J) 

AM  <2 . I  l*Y< I+lO.JI 
2  AM  <3*1  l*Y< 1+1 3. J> 

CALL  MTRXL  CAM.F.l.E) 

OYI 17. J)»Y<20. J) 

0YC18. J)*Y<21  .J) 

OYI 19. J)»Y<22. J> 

CALL  MTRXL<AM,0Y<1 7. J) . 1 . Y<2. Jl ) 
UP«Y<2. JI-WX 
VP*Y<3. JI-WY 
HP. Y<4. J1-9Z 

VA-SORT  <UP**2+VP**2+WP**2> 

SB-VP/VA 

SA«WP/VA 

QP*  < .5*RH0*VA  >*VA 

RM*VA/VS 

CALL  TABL < 3. J I 


0S730 

05740 

05750 

05760 

05765 

05770 

05760 

05790 

05800 

05810 

05820 

05830 

05840 

05850 

05860 

05870 

05880 

05890 

05900 

05910 

05920 

05930 

05940 

05950 

05960 

05970 

05980 

05990 

06000 


06010 

06020 

06030 

06040 

06050 


06090 
06100 
061 10 
06120 
06130 
06140 


S3 


n  n 


CALI-  MTRXL  < AM. Y (5. J ) . 1 . WDOT ( 1  )  > 

DV*REFL*.5/VA 

PO»P*DV 

CD*CNP»PO 

D0-PD»SB 

CY=>  +  CNA*SB-CD*SA 
CZ<=+CNA*SA  +  CD*SB 
GO  T0<200. 200.2  1.200. 200). J 

200  CALL  T ABL  <  4 . J 1 
C 

C  COMPUTE  TOTAL  AERODYNAMIC  FORCE  AND  MOMENT  COEFFICIENTS 
C 

C  MAL*1  INPUT  CENTER  OF  PRESSURE  .  MAL=?  'NPUT  RESTORING  MOMENT  COEFrIC 
C 

20 1  GO  TO  (3.4  )  .MAL 
3  CP«CMA 

CMA»-CNA*A8S ( CG-CP ) /REFl 
4  CMY»CMA»SA+CMQ* 3*DV 

CMX«CLO*F  1  NC-f  ( --  »REFL*  .5/VA  )  »CLP 
CMZ«CMQ#R*OV-CMA*SB 
CT*TSL+AE* (PSL-PA ) 

COMPUTE  ACCELERATIONS  IN  BODY  SYSTEM 
C 

IF  (KAN.E0.2I  CT«0. 

C0S2L  =  (Y(17.J)*Y(17.J)+Y(18.J)*Y(ia.J)-Y(19.U)*Y(19.J)  )/ <HR0**2 ) 
G=GCON1/HRO**3* ( 1 . -3 • *REQ*REU*GNU* < 1 .-3. *C0S2L )/HR0**2 ) 

CALL  MTRXL  ( AM . Y ( 1 7 . J ) . 1 . F ) 

FQA»QP*DSQ 

fqr-foa*refl 

FZ*F  Y 

0 Y ( 2 . J ) =  -F ( 1 1*G+ <CT-CX*F uA >/FM 

DY ( 3 . J ) =  -F (2 )*G- (CY*FQA)/FM 

DY(4.J>»  — F ( 3 ) *G- ( CZ*FQA )/FM 

WDOT ( 4 >  =  < (FY-FZ)#Q*R+CMX*FQR)/rX 
(HOOT  (5  )  *  (  (FZ-FX  )  #P*R  +  CMY*FQR  l/FY 
HOOT (6 ) = ( (FX-FY )*P*G+CMZ*FOR I/FZ 
1 F ( KAN. EQ . 2  I  GO  TO  41 
WDOT (4  ) =wDOT ( 4 I-DMKOMK (  t  )  *P 
WDOT  (  5  )  =  WDQT  (  5  )  -UMKOMK  (  2  )  *Q 
WDOT  (6  )  =  WD0T  <6  )  DMKOMK  t  2  )  *R 
50 1  DAMP=FL0W*(XLT-CG>**2 

WDOT  <  5 ) =  WD0T ( 5  1  +DAMP*Q/FY 
WDOT (6 ! «WD0T (6  »+DAMP*R/fZ 
41  CALL  MTRXL ( Y(8. J ) .WDOT (4 ) . 1 .DY (5. J) » 

DO  5  KA=8. 10 

KD»KA+3 

KE»KA+6 

DY (KA . J )=R*Y (KD. J  >~Q*Y (KE. J ) 

DY(KD. J)-P*Y<KE.J)-R*Y<KA. J) 

5  DY (KE . J >=Q#Y (KA, J >-P*Y (KD, J > 

CALL  MTRXL  < Y < 8 . J ) , DY . 2 , J ) , 1 , DY ( 20 . J  )) 

RETURN 


06150 
06160 
06170 
06180 
06190 
06200 
06210 
06220 
06230 
06240 
06250 
06260 
06270 
06230 
06290 
06300 
0631  0 

06330 

06340 

06350 

06360 

06370 

06380 

06390 

06400 

06410 

06450 

06460 

06470 

06480 

06490 

J6500 


06540 

06550 

06560 

06570 


06580 

06590 

06600 

06610 

06620 

06630 

06700 
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non  o  o  r>  o  n  n  n  nono 


END 


06710 


SIE  IT 

$ I  BET  C  MTR  LIST. REF. DECK 

SUBROUT  I NE  MTRXL ( A 1  , A2 . NC2 . °R0 ) 

DIMENSION  A1I3.3I.A2 (3. 3> .PRO (3. 3) 

DO  1  I  =  1  *  3 
DO  1  J=1,NC2 
PRO ( I , J )*  0.0 
DO  1  K= 1 . 3 

1  PRO(I.J)  -  PRO< I .Jl+Al ( I ,K)*A2<K.J) 

RETURN 

END 

SIBFTC  CHK  L 1  ST .REF, DECK 

SUBROUTINE  CHKOUT 

CC(-  N  Y <22.5 ) ,DY<22.5 ) .R0FF1 (22.5 ) 

IN  ALL  COMMON  STATEMENTS  D=DR I VER . A= AOPTUN « B=BFACT . E=EVAL . C=CHK , I = 1 NT 
R2=RKERR2.R=RKG.T=TABL.F=PhASIN,P=P1NTRJ.X=X0UT. AND  L=T I2LP 

COMMON/DCP/  TO 

OMMON/OC/  NUF.IOUT 

COMMON/DCF /P 1 . TBO 

COMMON/DC  1/  H 

COMMON/EC./  CHEK.RS  , 

COMMON/CF/  TFP 
COMMON/PC/  ZIM 
COMMON/DEC/EPSQ 


1 OUT -  1 ,  NO  OUTPUT — I  OUT  =2 .  PRINT  OUTPUT 
C  1  OUT  *  3 .  END  OF  PHASE — I0UT=4«  IMPACT 
C 

NUFrO 

lFIYll ,5) .EQ.TFPIGO  TO  5 

CHEK  =  S0RT  < Y  < 1 7.5 )**2+Y< 18,5  )**2+Y ( 19,5  1**2 ) 

RS=20926428./ < ( 1 .+EPSQ*<  Y< 1 9.5 l/CHEK )**2/ < 1 .-EPSQ) >**.5 1 
EMPT  *RS+Z I M 
IFICHEK.LE.EMPT )NUF»1 

IFINUF.EQ. I .AND.YI 1 ,51 .GT.20. 1G0  TO  2 
ENDPH= ABS ( TBO— Y (1,51) 

IFIENDPH.LT. .001 )G0  TO  3 
D I F*  ABS  <  Y ( 1  ,5)-T0) 

PO*AMOD<DIF.Pl 1 

IF<ABS<PO-PI 1 .LT. .0001 1GO  TO  4 
IFCPO.LT..001 )G0  TO  4 
PO*P!-PO 


06720 

06730 

06740 

06750 

06760 

06770 

067B0 

06790 

06800 

06810 

06820 

06830 

06840 

06850 

06860 

06870 

06880 

06890 

06900 

06910 

06920 

06930 

06940 

06950 

06960 

06970 

06980 

06990 

07000 

07010 

07020 

07030 

07040 

07050 

07060 

07070 

07080 

07090 

07100 

07110 

07120 

07130 

07140 

07150 

07160 

07170 

07180 

07190 

07200 

07210 

07220 


SS 


H- AM INI (H.PO.ENDPH) 

07230 

! OUT  « t 

07240 

GO  TO  1 

07250 

4  I  OUT-2 

07260 

H-AMINi (H.Pt .ENDPH) 

07270 

GO  TO  1 

07280 

3  I  OUT-3 

07290 

T  0«  TBO 

07300 

GO  TO  I 

07310 

2  ! OUT-4 

07320 

GO  TO  1 

07330 

5  TO-TFP 

07340 

endph-tbo-tfp 

07350 

GO  TO  4 

07360 

l  RETURN 

07370 

END 

u  7380 

*IBFTC  INTEGR  L 1ST .REF . DECK 

07390 

SUBROUTINE  INTEG 

07400 

COMMON  Y < 22 «  5 )  « DY ( 22 • 5 ) • ROFFI  <22«5) 

07410 

C 

07420 

C 

IN  ALL  COMMON  STATEMENTS  D-DR I VER . A-AOPTUN, B-BFACT . E-EVAL . C-CHK . 

5-  INT 

0  7430 

c 

R2-RKERR2.R=RKG.T=TABL.F-PhAS1N.P=PINTRJ,X*X0uT.AND  L=TI2LP 

07440 

c 

07450 

COMMON/DC  I /  H 

07460 

c 

07470 

COMMON/DI/  DONE 

07480 

c 

07490 

COMMON/D I X/  TIME 

07500 

c 

07510 

COMMON/D  I R2/  NOGOOO 

07520 

c 

07530 

COMMON/B I R2/  YSAVE ( 22 ) . QYS A VE ( 22  > • ROFFS ( 22 ) 

07540 

c 

07550 

COMMON/ I R2/  IND. YHALFI22 ) 

07560 

c 

07570 

c 

INTEG  IS  THE  INTEGRATION  MONITOR  WHIC  CONTROLS  THE  RUNGA  KUTT  INTEGRAT 

07580 

c 

AND  THE  ERROR  CHECK 

07590 

c 

SAVE  THE  INITIAL  CONDITIONS 

07600 

DO  I  1=1.22 

07610 

YSAVE ( 1 )*V< 1 .51 

07620 

DYSAVE  < I  >»DY( I . 5  > 

07630 

I  ROFFS ( I > -ROFFI (1.5) 

07640 

CALL  RKG(H.TIME) 

07650 

3  NOGOOD-O 

07660 

CALL  RKERR2 ( H ) 

07670 

I F ( IND.EQ. 19)00  TO  24 

07680 

I F ( NOGOOD • EQ • 0 )  GO  TO  23 

07690 

c 

IF  NOGOOO  EQUALS  0  THE  INTEGRATION  IS  ACCEPTABLE.  IF  IN  ADDITION 

IND  I 

07700 

c 

EQUAL  TO  19  DT  IS  DOUBLED  FOR  THE  NEXT  PASS. 

07710 

IF (H.LT. .0001 1  GO  TO  25 

07720 

H-H/2. 

07730 

Y(|  .1  )*TI ME 

07740 

DO  2  1=2,22 

07750 

56 


2  Y( I ,5)=YHALF< I ) 

CO  TO  3 

25  DONE*2.0 

:  DONE  IS  AN  INDICATOR  TO  TELL  IF  THE  SIMULATION  SHOULD  BE  CONTINUED  l.Y 
RETURN 
24  H*2.0*H 
23  DONE* 1 .0 
RETURN 
END 

> IBFTC  RKER2  L I  ST , REF . DECK 
SUBROUTINE  RKERR2  <H> 

COMMON  Y  <  22 . 5 ) •  DY ( 22 * 5 ) , ROFF 1  <  22 . 5  > 

IN  ALL  COMMON  STATEMENTS  0*DR I VER . A= AOPTUN . B*BFACT , E=EVAL , C*CHK * I  * INT 
R2=RKERR2*R*RKG,T=TABL.F*PHASIN,P=PINTRJ*X=XOuT.AND  L=TI2LP 

COMMON/B I R2/  YSAVE ( 22 ) *DYSAVE  <22 ) ,ROFFS<22 ) 

COMMON/ I R2/  IND. YHALF <22 ) 

COMMON/FR2/EPT I  NY . EPB I G 

COMMON/D  I R2/  NOGOOD 

DIMENSION  Y  I  (22) *TRUNC<22 >  * YI I (22) 

Y I (1  ) *Y  < 1 ,5) 

DO  1  1=2*22 
YI (I )*Y< I ,5) 

1  CONTINUE 

DO  7  K=2  *  22 
Y  <K  *3  > = YSAVE  <K ) 

DY(K.5)*DYSAVE<K) 

ROFFI <K.5)»R0FFS(K) 

7  CONTINUE 
I  ND*  1 
HALF*H/2. 

CALL  RKG  (HALF* Y< l . I  )  ) 

DO  10  K*1  ,22 
10  YHALF  <K)*Y(K*5) 

CALL  RKG  (HALF* Y( 1 *3) I 
YI I (1 >=Y< 1 .5) 

DO  3  I  1*2,22 
YI  I  ( I  I  )*Y( I  I ,5) 

3  CONTINUE 

DO  4  1=5.22 

IF(ABS< YI  I < I  )-YI  < l  )  > *LT. .000001  )  GO  TO  2 
TRUNC ( I  1**02222222  *ABS < (YII(I)-Yl(I) >/Y 1(1)1 
I F  <  TRUNC ( 1  ) *  GT  *  EPB I G I GO  TO  5 
IF  (TRUNC  (  I  )  .LT.EPTINY  )  INOMND+I 
GO  TO  4 

2  INOMND+l 

4  CONTINUE 
GO  TO  6 


07760 

07770 

07780 

07790 

07800 

07610 

07820 

07830 

07840 

07850 

07860 

07870 

07880 

07890 

07900 

07910 

07920 

07930 

07940 

07950 

07960 

07970 

07980 

07990 

08000 

08010 

08020 

08030 

08040 

08050 

08060 

08070 

08080 

06090 

08100 

08110 

08120 

08130 

08140 

08150 

06160 

0817Q 

08180 

08190 

08210 

08220 

0823Q 

08240 

08250 

08260 

08270 

08280 


1 - 
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5  N0G00D=1  08290 

6  RETURN  08300 

END  0B310 

4IBFTC  RKGA  LIST* REF • DECK  08320 

SUBROUTINE  RKG  (H.TIME)  08330 

COMMON  Y(22»5>.DY(22.5>,H0FF1 (22.5!  08340 

DIMENSION  A <4 ) .8(4 ) ,C <4 )  08350 

DATA  All ).B<1 >*C(1 ).A(2).8(2).C*2).A(3).B(3>.C(3).A(4).B(4).C(4>/  08360 

1 .5.2. .  .5. .2928932. 1 . . .2928932. 1 . 7071068.1 . . 1 .7071068. . 1 6666666.2.  08370 

2..  5/  08380 

J>  i  08390 

99  GO  TO  ( 1 .2.3.4) . J  08400 

1  Y ( l . 1  ) =T I  ME  08410 

DO  11  1=2.22  08420 

Y(  1 . 1  )*Y (  I  ,5)  08430 

0Y< I . 1 >*DY< I .51  08440 

ROFF 1(1.1)  =ROFF 1(1.5)  08450 

II  CONTINUE  08460 

GO  TO  5  08470 

2  Y < 1  ,2 ) =Y ( 1  .  S  )+H/2.  08480 

GO  TO  6  08490 

3  Y( 1 ,3)=Y( 1 .2)  08500 

GO  TO  6  08510 

4  Y( 1 ,4 )=Y( 1 ,3>+H/2,  08520 

6  CALL  EVAL  <J>  08530 

5  DO  50  1=5.22 

Y ( I , J+l )=Y( l . J)+H* (A( J>* (DY ( I . J)-B< J )*«0FF1  (I.J>)>  08550 

ROFF 1 < 1 » J+l >  *  R0FF1 ( I . J )+3»* ( A ( J >* ( DY ( I . J ) -8 < J > *ROFF 1 ( I . J ) > )  08560 

1  — C ( J ) *DY (  1  . J )  08570 

50  CONTINUE  08580 

l F ( J • EQ . 4 ) GO  TO  52  08590 

GO  TO  53  08600 

52  JT*JFl  08610 

Y ( 1 ,5 ) =Y ( 1 .4 )  08620 

CALL  EVAl(JT)  08630 

GO  TO  100  08640 

53  J»JFl  08650 

GO  TO  99  08660 

100  RETURN  08670 

END  08680 

4I8FTC  TAB  I  ST. REF. DECK  08690 

SUBROUTINE  TABL (LOOKUP. J )  08700 

COMMON  Y < 22 .5 ) . DY (22 . 5 ) .R0FF1  < 22 . 5  )  08710 

c  08720 

C  IN  ALL  COMMON  STATEMENTS  D  =  DR I VER . A  =  AOPTUN. B  =  BFACT . E =E V AL . C*CHK . I  *  I  NT  08730 

C  R2»RkERR2.R=RKG.T=TABL.F=PhASIN.P=PINTWJ.X*X0UT,AND  L=TI2LP  08740 

c  08750 

COMMON/PT/  N.HHI44 ) ,RH0<44 ) , VS (44 ) .PA (44 )  08760 

c  08770 

COMMON/OAETX/  WXS.WYS  08780 

c  08790 

C0MM0N/D8TP/NM . WP ( 50 )  08800 

c  08810 

S8 


oo  oooooo  no 


COMMON/DT/  11.12.13.14  08820 

C  08830 

COMMON/ ABE TP/  WXP ( 50 » . WYP ( 50  I  08840 

C  0B'.'50 

COMMON/BET/  HT  08860 

COMMON  /FT/  PM(20 i .CX (20 ) .CNA (20 ) .CMA (20 ) .CMGi(20 ) . T (40 ) . FX (40 ) .  08870 

1FY (40 ) . TSL (40 > . CG ( 40 ) . FM ( 40 i • NM . NT . CND ( 20 ) .CNP(2Q ) 

C  08880 

COMMON/ET /DENS iSPD. PRES . CD . CL . CP . CMT . TME ( 5 ) .RM.DMKOMK (3 ) .FLOW  088 

1 .CLD.CLP 

08900 
08920 

COMMON/EFT/KAN  08930 

08940 

I  1.12.  13.  14  IS  WHAT  LAST  ENTRY  FROM  THE  ATMOS . W 1 ND . MACH  AND  TIME  TAB  08950 

08960 

GO  TO  ( 1 .2.3.4) .LOOKUP  OP  0 

08980 

INTERPRET  FOR  ATMOSPHERIC  PARAMETERS  08990 

09000 

1  1 F ( HH (Il).LT. HT )  GO  TO  5  09010 

7  IF(HH<  I  1-1  ) .LT.HT )  GO  TO  6  09020 

IF  (I1.EQ.2)  GO  TO  6  09030 

11*11-1  09040 

GO  TO  7  09050 

5  IF  (Il.EQ.N)  GO  TO  6  09060 

11*11+1  09070 

IF  (Il.LT.N)  GO  TO  1  09080 

6  CON* ( HH (11-1) — HT ) / ( HH (  I  1  - l  ) -HH (II))  09090 

DENS  =  RHO(  I  1-1  )  —  ( RHO (  I  1 - 1  )-RHO(  I  1  )  )*CON  09  100 

SPO  =  VS( 11-1  )-<VS( I  1-1  >-VS<  I  1  )  >*CON  09110 

PRES*PA( 1  1-1  )- (PA ( I  1-1  ) —PA (  1  1  )  >»CON  09120 

RETURN  09130 

09140 

INTERPRET  FOR  W  I  NO  09150 

C  09160 

2  IF  (NW.EO.l )  TO  TO  15  09170 

14  IF  (HP( 12). LT.HT)  GO  TO  10  09180 

12  IF  <HP( 12-1 ). LT.HT >  GO  TO  11  09190 

12*12-1  09200 

IF  (I2.EQ.1)  GO  TO  17  09210 

GO  TO  12  09220 

10  IF  (I2.EQ.NW)  GO  TO  16  09230 

12*12+1  09240 

GO  TO  14  09250 

11  CON* (HP( 12-1 >-HT >/<HP< 12-1 )-HP( 12 ) )  09260 

WXS»WXP< 12-1 >- (WXP( 12-1 )-W.<P( 12 ) >*CON  09270 

WYS*WYP< 12-1 ) — ( WYP ( 12-1 ) -WYP ( 12) )*CON  09280 

GO  TO  13  09290 

17  12*2  09295 

16  WXS*0,  09300 

W YS*0 .  09310 

GO  TO  13  09320 
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IS  KXS*KXPUI  09330 

KYS-WYP(l)  09340 

1 3  RETURN  09350 

09360 

INTERPRET  FOR  PARAMETERS  VERSUS  MACH  09370 

09380 

3  IF  (PMI I3I.LT. RMIGO  TO  20  09390 

22  IF  (PMI 13-1 J.LT.RM)  GO  TO  21  09400 

IF  (I3.EQ.2)  GO  TO  21  09410 

I 3» 1 3- 1  09420 

GO  TO  22  09430 

20  IF  ( I 3.EQ.NM )  GO  TO  21  0944Q 

13* 1341  09450 

! F  ( I 3.LT.NM )GO  TO  3  09460 

21  CON*(PM( 13-1 )-RM)/(PM( 13-1 )-PM( 13)  I  09470 

CD»CX( 13-1 >-<CX< 13-1  )-CX  < 13) )*CON  09480 

CL»CNA< 13-1  )-<CNA( I  3—1  )-CNA( 13) )»CON  09490 

CP*CMA< 13-1  )  —  ( CMA  < 13-1  )~CMA(I3) )*CON  09500 

CMT*CMQ( 13-1 ) — ( CMQ  < 13-1 )-CMQ<  13)  >*CON  09510 

CLD»CND ( 13-1 )  —  ( CND  < 13-1  )-CND< 13) >*CON 

CLP*CNP( 13-1 ) — ( CNP ( 13-1 )-CNP, 13) )*C0N 

RETURN  09520 

09530 

INTERPRET  FOR  PARAMETERS  VERSUS  TIME  09540 

09550 

4  IF  (NT.EQ.l)  GO  TO  40  09560 

32  IF  <T  I4I.LT. Y(l.JI)  GO  TO  30  0957o 

IF  <T(  14-1  I.LT.YIl .J)  )  GO  TO  31  09580 

IF  (I4.EQ.2)  GO  TO  31  0959C 

I4« 14-1  09600 

GO  TO  32  09610 

30  IF  (I4.EQ.NT)  GO  TO  33  09620 

14 ■ 144 |  09630 

IF  (  I4.LT.NT  )C,0  TO  32  09640 

31  CON* <T< 14-1 >-Y( 1 ,J> >/(T< 14-1 )-T < 14 ) )  09650 

TMEI  1  >*FX< 14-1  )-<FX( 14-1  l-FX(  14 ) )  *CON  09660 

TME  ( 2  )  *FY  (  14-1  )-(FY<  14-1  )-FY(  14  )  )*CON  096"’0 

TME(3)*TSL( 14-1  1  ( TSL ( I  4- 1  ) - TSL < 1 4 ) ) *CON  09o80 

TME (4  ) *CG ( 14-1  )- (CGI  14-1  )-CG(  14)  )*CON  C9690 

TME ( 5  ) *FM ( 14-1  )-(FM( 14-1  )-FM(  14 )  )»CON  09700 

34  DMXOMK ( 1  ) a ( FX ( 1 4  1  ) -FX ( 1 4 )  )/ ( T ( 1 4- 1  ) -T ( 1 4  )  ) /FX  09720 

OMKOMK ( 2 ) ■ ( F Y (  I  4  l  > -F  Y  (  |  4  )  >/ (  T  (  1 4- 1  ) -T  (  1 4  )  )  ,'FY  09730 

OMKOMK (3 ) aDMKOMK (2 )  09740 

FLOW- <FM(  14-1  ) — FM (  14)  ) / ( T (  14—1  )  —  T<  14)  ) 

GO  TO  13  09750 

33  TME ( 1 ) «FX (NT )  09760 

TME(2 >*FY(NT )  09770 

T  ME ( 3 ) *  T  SL ( NT )  09780 

TME (4 ) *CG ( NT )  09790 

TME (5 )*FM(NT )  09800 

IF  (KAN. £0.2 )  GO  TO  201  09820 

GO  TO  34  09830 

40  TME  (  1  l-FXU  )  09840 
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n  o 


i 
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TME(2)=FY(1)  09850 

TME ( 3 ) =  TSL ( 1 >  09860 

TME  ( *  )  =CG  ( 1  I  C'9870 

TME(5)=FM(1)  09880 

201  DO  200  1*1.3  09900 

200  DMKOMK (I ) =0 .  09910 

RETURN  09920 

END  09930 

* IBFTC  FAZ1N  L 1  ST . REF . DECK  099*0 

SUEROUTINE  PHASIN(J)  09950 

09960 

IN  ALL  COMMON  STATEMENTS  D»DR I VEK . A * AOPTUN . B=UFACT . E*E VAL . C  =  CHK . 1  *  1 NT  09970 
C  R2*RKERR2.R*RKG.T*TA8L.F=PhASIN.P»PINTRJ.X*X0UT.AND  L=T I2LP  09980 

C  09990 

COMMON/DCF/P 1 .TBO  10000 

C  10010 

COMMON/DF/  TMII  10020 

C  10030 

COMMON/CF/  TFP  100*0 

C  10050 

C0MM0N/FR2/EPTINV.EPB1G  10060 

C  10070 

COMMON  /FT  PM ( 20 ) • CX  <  20 ) . CNA  <  20 ) *  CM A ( 20 ) • CMQ (20)*T (40).FX(40) .  10080 

1FY (40 ) «  T  SL  <  40 ) • CG ( 40 ) . FM ( 40  ) .  NM.NT.CND(20 > .CNP(20 ) 

C  10100 

COMMON  /PF/  TBOB( 10 ) ,AEB< 10 ) .OSQB( 10) .REFLBl 10) . TMI IB( 10 > .  10110 

1 P I  8 ( 1 0 ) . EPT I N0 ( l  ) .EPBIGBC 1 0) .TFP8( 1 0) .KANB( 1 0) .  NMB<10)  10120 

2 . P MB <  20 , 1 0 ) , CXB ( 20 . 1 0 ) . CNA8 (20.10). CMAB (20.10). CMQB (20.10).  10130 

3NTB ( 1 0 ) . TB(40 . 1 0 ) .FX8(4Q . 10 ) .FYB(40. 10). TSLB(40. 10 ) .Cl»o.40. 10).  101*0 

4F MB  (40*20  )  .XLT  (  1  )  ,CNDB(20. 10)  .CNPB(20.  10  )  .FINCBUO  ) 

C  10160 

COMMON/EFt/KAN  10170 

C  10180 

COMMON/EF/  REFL . AE . DSQ. XLNTM.F I NC 

C  10200 

TBC»TBOB(J>  10210 

AE«AEB(J)  10220 

DSO»DSOB(J)  10230 

REFL*REFLP ( J)  102*0 

KAN«KANB(J)  10250 

TMI I-TMI IB( J)  10260 

P|*PIB<J)  10270 

EP8 I G«EPB I G8( J )  10280 

EPTINY*EPTIN8< J)  10290 

F 1 NC*F I NCB  <  J ) 

TFP-TEPB(J)  10300 

XLNTM*XLT(J>  10315 

NM«NM8(J>  10320 

DO  200  1 ■ 1 .NM  10330 

PM< | )«PM8( I . J)  103*0 

CX( I  )»CXB( I • J)  10350 

CNA(  1  )»CNAB(  I  *  J)  1  03*>0 

CM A ( 1  )  »CMAB ( | . J )  10370 
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CMO( I >»CMOB( I . J)  10380 

CNB ( I  l*CN06 ( 1 . Jl 

200  CONTINUE  •  i 0390 

IF  (KAN.E0.2)  CO  TO  202  ;0400 

201  NT^NTBIJ)  10410 

00  203  I  a  I .NT  10420 

T ( I )«T8( I . J>  10430 

FX< I l*FXB( I.Jl  10440 

FY( I )«FYB< 1 • Jt  10450 

T  SL  < I  I'TSLBl I «J>  10460 

CG ( 1  I *CGB <  I • J I  10470 

FM( I )*FMB( I . Jl  10480 

203  CONTINUE  10500 

202  RETURN  10510 

END  10520 

3I0FTC  P1NTR  L I  ST .REF. DECK  10530 

SUBROUTINE  PINTRJ  10540 

C  10550 

C  IN  ALL  COMMON  STATEMENTS  0»0R 1 VER. A* AOPTUN. JaBFACT . £»EV AL . C»CHK . I  =  I  NY  10560 

C  R2»RKERR2.R=RK0.T*TABL.F*PmAS1N.P=PINTRj,x*X0UT.AND  L-TI2LP  10570 

C  1 0580 

COHXON/PT/  N.HH ( 44  I • RHO < 44 )  . VS (44 ) . PA (44  I  10590 

C  10600 

COXMON/DAP/TM . AL. XMANT . YKANT . DELC • DELT .DELH.TOWT I L  10610 

C  10620 

COMMON/CAEP/NTYPE  10630 

C  10640 

COMMON/DBTF/NW.HP.50 >  10650 

C  10660 

COMMONZDCPZ  TO  10670 

C  1 0680 

COXHOn/AP/  XF IRST . VF IRST • XCHANG . YCHANG . XL AST . YLAST . NBF . 0F (501.  1 0690 

10ELTA  10700 

C  10710 

COMMON/DP/  UO.VO.KO.ZL.XLAT .NFAZE.NBST. IFOUT  10720 

C  10730 

COMMON/ABETP/  KXP100 ) . KYP(50)  10740 

C  10750 

COMMONZEPZ  MAL  10760 

C  10770 

COMMON/’PCJ  ZIM  10780 

C  10790 

COMMON/BPX/  I PEN  10800 

C  10810 

COMMON  /PF/  T60B< 10  I , AEB< 10  I .DSQB ( 1 O) .REFLBI 101 . TM I  IB< 10 ) .  10820 

1 P I B ( 1 0  I . EPT INB  < 1  I «EP0 I GB( 10  1. TFPB (10) .KANB (101.  NMB(IO)  10830 

2 . PMB (20.10)*CXB(20 • I  0 1 . CNAB(2C .10). CMAB ( 20 . 10 1 . CMOS (20.10).  10840 

3NT0 (10>. TB( 40. 10). FXB (40.10). FYB (40.10) .TSLB(40. 10 ) .CG8(40. 10 ) .  10850 

4FMB (40.20 ) . XL T ( 1  ) .CNDd<  20 . 10 ) .CNPB (20 .10) .F I NCB ( 10  > 

C  10870 

01  MENS  1  ON  DMACH (6 ) .DTI  ME  <  5 1 . 1  DEN ( I  i , CND<20 ) .CNP(20 ) 

dimension  PM (20 ) ,CX (20 ) .CM'  ,2,1  • .CM A (20 ) .CMQ (20 ) .T (40  > .  10890 

1FX(40).FV(40)  .TSL( 40)  .CG( 40  ).F  -  1 


non  non  nor.  n  n 


10910 

NT  VPE  -TYPE  OF  OUTPUT  DESIRED  <  *  1  . REGULAR  TRA J > . < *2 . I NTERATES >  10920 

(=3. ANGULAR  OUTPUT ). <«4. BALL  FACTORS >.( »5 .PARAMETER  VAR  1 AB I L 1 T Y .D )  10930 

I FOUT a 1  CONTINUE  PRINTING  AT  END  OF  PHASE-NO  SKIP  TO  NEXT  PAGE  10940 

10990 

1  READ  <5.1  COO ) NT/PE . 1 FOUT  10960 

RE AD <5. 1500)1  DEN  10970 

WRITE (6. 499) I  DEN  10980 

10990 

NFAZE  IS  THE  TOTAL  OF  THE  PHASES-  NBST  IS  PHASE  TO  PICK  BOOSTER  UP  AT  11000 

11010 

READ  <5.2000)  T0.2L.2IM.XLAT.XL0NG.NFAZE.NBST  11020 

WRITE  (6.900)  T0.ZL.ZIM.XLAT.XL0NG.NFA2E.NBST  11030 

J=1  11040 

11050 

MAL*I  INPUT  CP  MAL*2  INPUT  CMA  11060 

I  1070 

READ(5.7000)UO.VO.WO.MAL  11080 

WRITE<6.510)UO.VO.WO.MAl  11090 

read<s.3loo)th.al.xwant,ywant.oelta.towtil  i i ioo 

WRITE<6.S20)TH.AL.XWANT,YWANT.0£LTA.T0WT!L  1 1 1  1C 

READ <5.2500 )N. <HH< I ) .RHO< I ) . VS < 1 ) . PA < I > . I » 1 . N >  11120 

WRITE <6. 530 )N , <HH< I  ) .RHO. I  ) • VS ( I  ) .PA  <1  ) • I  *  1 «N )  11130 

I F (NTYPE .NE.3 )G0  TO  2  11140 

READ<5.3900 )XFIRST.XLAST.XCHANG. YF 1 RST . YLAST . YCHANG  11150 

WRITE<6.540)XFIRST.XLAST.XCHANG.YFIRST.YLAST.YCHANG  1  1  160 

GO  TO  1 2  11170 

2  KEAO<5.3500 >NW. < HP < I ) .WXP< I ).WYP< 1 ) . I »1 .NW )  11180 

WRITE  <6.550  >NW.  < HP < I  ) . WxP< I  ) . WYP < I  ) •  I >1 .NW  >  11190 

IF(NTYPE.NE.2)G0  TO  4  11200 

RE AD <5. 4 000 )NBF. <BF< I ) . I ■ 1 . N8F )  11210 

WRITE  <6.560 )N8F,  <0F ( I  )«  |s]  . NBF )  1  1220 

12  READ  <  5 . 3000 )DELC .CELT . DELH  11230 

WRITE  <6.570 )DELC.OELT.OELH  1 1240 

c  11 250 

C  TFP.T I  ME  OF  FIRST  PRINT  KAN* 1 .MACH  AND  TIME  TABLES  KAN*2 « MACH  ONLY  11260 

C  I OPTUN«  30  TO  GROUND  11270 

c  1 1 280 


4  READ  (5,4500)  TBO, AE  »DSQ .REFL . XLENTH, EPT I  NY ,EPB I G. F INC .TMII.PI.TFP 
I  •  KAN 

WRITE  <6.580)  TBO.AE.DSO.REFL.XLENTH.EPTINY.EPBIG.FINC.TMI  I .PI .TFP 
1  ,  KAN 

READ  (5.5000)  NM . <PM ( I ) • CX  < I ) ,CNA  < I  ) ,CMA< I  ) .CMOI  I  >,CND< I  I . CNP < I  ) , | 
1 = 1 .NM ) 

WRITE  <6.590 )  NM. (PM< I ) ,CX< I ) . CNA ( I ) .CMA ( I ) . CMQ< I ) , CND < I ) ,CNP< I ) , I 


l » 1 .NM | 

I F (NT YPE .NE .5  >G0  TO  7  11330 

RE AO  <  5  «  5500 ) <  OMACH  < I  ) « I  a  1 , 6 ) 

WRI TE <6.600 ) (OMACH < I  ) • I  *  1 < 6  I 

DO  8  K=1 ,NM  11360 

CX<K >-CX<K)*DMACH< 1 )  11370 

CNA (K ) *CNA  <K ) *DMACH  <  2  >  11380 

CMA(K>*CMA(K)*DMACH(3)  11390 


(S3 


CMG  IK ) *CMQ 1 K ) *0MACH 1 4 ) 
CND ( K  >  »CND ( K  > *OMACH ( 5  I 
CNP  <  K  )  «CNP ( K ) »OM ACH <  6  > 


10 


15 

9 


20 

6 


499 

500 
510 
520 


T808I J)«T80 

11410 

AE8< J  t*A£ 

1  1420 

OSQB1 J)*OSQ 

1  1430 

REFLB! J)*REFL 

11440 

TM| 1R( J)*TMI  1 

1  1450 

P1B< Jl-Pl 

1  1460 

EPT1NB1 J)*EPT1NY 

1  1470 

EP8IGB*J)«eP8IG 

1  1480 

FINCBt  J)*F1NC 

TFPB! JI-TFP 

11490 

KAN8  (.)>«*  AN 

1  1500 

alt t j»«xlenth 

11515 

NMB 1 J ) *NM 

11520 

DO  10  1«1 .NM 

1  1530 

PMB1 1 «J)*PM< | 1 

11540 

CXBt T . J)*CX1 I » 

1  1550 

CNAB( 1 . Jl-CNA! 1 ) 

11560 

CMAB< 1 . J)«CMAI | 1 

1  1570 

CXOBI 1 .  Jl-CMOl 1 1 

CN081 1 . JI«CMQ( 1 1 

CNPBl 1 . JI*CNPl I > 

IFIKAN.NE.I )GO  TO  6 

1  1590 

BE  AO  1  5.  6000  >NT.  1T1  1  ).FXU  )«FY1  I  »  <T5L(  1  I.CGU)  .EMI  1  >  . 

1*1. NT ) 

1  1600 

WRI TEl 6.610 )NT. ( T ( 1) .EX ( I  ) , FY ( 1  1.TSL1 1  ) . CG < 1  ) .FM 1 1  ) , 

1*1 .NT ) 

1  1610 

IFlNTYPE.NE.5lGO  TO  9 

1  1620 

RE AD (3. 65001 (DT 1 BE  *  1  » . I *1 .5  1 

WRITE (6. 620) 1DTIME1 I ) . 1 «l .5 1 

DO  15  K*1 .NT 

1  1650 

FX IK  >»FX 1K)*DT 1ME! 1  t 

1  1660 

FY IK )*FY (K )*0T IHE12 ) 

11670 

TSL1K)*TSL<K)*0TIME.  13) 

11680 

CG  (K 1  =  CG 1 K ) *DT | ME  1 4 ) 

11690 

FM1K)*FM1K)*0TIME<5> 

11700 

CONTINUE 

NT01 J)=NT 

11720 

DO  20  1=1 .NT 

11730 

Tail . J)»T 1 1 ) 

11740 

FXB1 1 .JI«FX( I ) 

11750 

FYBt I . J)«FY1 I ) 

1  1760 

TSLBt I . J ) =  T  SL 1 1  ) 

1  1770 

CGBl 1 . J ) »CG( I ) 

1  1  780 

FMB 1 1 . J) *FM1 | ) 

1  1790 

CONT 1 NUE 

J«  J-M 

11810 

IF1J.LE.NFAZE 1 GO  TO  4 

1  1820 

RETURN 

l  1830 

FORMAT  1 IH1 • 10X. 1 2A6 ) 

1  1840 

FORMAT  1 1 IX. FI  1 .2.2F1 1 .2.2F1 1 .5.7X. 14. 7X. 14 > 

11850 

FORMAT! I 1X.3FI 1.  .7X.I4) 

1  1860 

FORMAT!  1  1X.4F1  1 .2. FI  1 «0«F1 i .1  ) 

11870 
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530  FORMAT ! 5X , 1 4 a 2X a F9.0 . 2X . F 1 1 .9 . 2F 1 1 . 3/ < 11 X . F9.0. 2X . F 1 1 . 9, ZF 11 . 3 )  )  1  I860 

540  FORMAT <1 1X.6FI 1 .  )  11  890 

550  FORMAT <5X, I4.2X.3F1 1 .2/0 SX.3F11 .2 J 1  11900 

560  FORMATtSX. I4.2X.F1 1 .6/! 1 1X.F1 1 .6) >  11910 

57G  FORMAT! 1  1X.3F1  1 .5)  11920 

580  FORMAT  ( 1  1 X  a F l  1 .2 , 2F l 1 . 7 . 4F 1 1 .4 .F 1  l . 9 . 2F 1  1 . 4/F 1 1 ,4 .  I  1  > 

5°0  FORMAT  (5X. I4.2X.F1 1 . 3. FI  1 .5. FI l .2. FI  1 .3*.  I  .2 , FI  1  . 3  .F 1  l  .  3/ (  1  1 X  ,F  l 
11.3iF11.5.FU,2.Fll,3.Fll.2.FU.3.F11.3)) 

600  FORMAT  ( 1  1 X  a  6F 1 1 .2  > 

610  FORMAT (5X a  I4.2X.F1  1 .2 .FI  1 .3, FI  1 . 1  .FI  1 .2a F 1 1 .4, FI  1 .3  119 

1  (llX.F11.2.Fll.3.Fll.l.Fil.2,FU.4.Fi;.3)) 

620  FORMAT! 1 1X.6F1 1 .3)  11990 

iOOO  F0RMAT(8X.2! II »7X1 1  12000 

1500  FORMAT ( 1 2 A6 )  12010 

2000  FORMAT  !8Xa5F8. 0.2! 12.6X1)  12020 

2500  F0RM_T<6X« I  2 a 4F8 . 0/ ! 8X , AF8 . 0 1  >  12030 

3000  FORMAT  <8X,6F8«Q)  1204© 

3500  FORMAT ( 6X .12. 3FS eO/ ( 8X, 3F8 * 0)1  12050 

4000  FORMAT !6X. 12.F8,  )  12060 

4500  FORMAT  < 1 0F8. 0/F8 . 0 ,  (11,7X1) 

5000  FORMAT  ( 6X a  1 2 . 7F8 . 0/ ( 8X , 7F8 . 0 )  ) 

5500  FORMAT  (16X.6F8.  ) 

6000  FORMAT (6X , I  2  a 6F6 . 0/ ( SX a 6F8 . 0 >  ) 

6500  FORMAT! 16Xa6F8aO)  12110 

7000  FORMAT (8X.3F8.0a 11 )  12120 

END  12130 

SIBFTC  XOUTP  ..1ST,  REF,  DECK  12140 

SUBROUTINE  XGUT  (Jr, IX)  12150 

COMMON  Y (22,5 ) . DY (22,5) ,R0FF1 (22.5)  12160 

C  12170 

C  IN  ALL  COMMON  STATEMENTS  D  =  OR 1 VER . A  =  AQPTUN, B=8FACT a E=EVAL . C  =  CHK a  1  =  INT  12180 
C  R2  =RKERRH  a  R  =  RKG  a  T  =  f  ABL  a  F  =  PHAS 1 N  a  P  =  P i NTR J  »  X  =  XOUT  a  AND  L  =  TI2LP  12190 

C  12200 

COMMON/XL/  TH,AH,D(3) .8(3! aV  12210 

C  12220 

COMMON/D AETX/  WXSaWYS  12230 

C  12240 

COMMON/D  IX/  TIME  12250 

C  12260 

COMMON/DX/  1JK  12270 

C  122B0 

COMMON/ABXL/  XYZ<3)  1229Q 

C  12300 

COMMON /BP X/  IOEN  12310 

C  12320 

DIMENSION  IOEN I  12) 

695  FORMA  T ( 23X  a  3F 15.2)  12330 

696  FORMAT (22X.3F15. 1 .2X.2Fl5.3l  12340 

697  FORMAT !/9XaF15.3aFl 3. 1 a2Fl 5. i .2X.3F 1 5.3)  12350 

698  FORMAT ( 1  HO , 1 7X  a  4HT I  ME  a l 2X a  1 HX . 1 3X a  1 Hy , 12X a . HZ . 1 7X.3HWX3. 10X.5HTHET  1 2360 

1  A.  1  IX.  5H  ALPHA /33X  a  2HXD.  12X.2H'i0a  1  1X,2HZD«  16Xa3HWYSa  1  OX  a  1HV/33X.3HX  12370 

200  a  1 1 X • 3H YDD a  1  OX  a  3HZDD )  12380 

699  FORMAT  ( ///40X40H 1 NTEGRAT I  ON  INTERVAL  IS  LESS  THAN  .00001)  12390 

6b 


no  no 


I F ( ! X.£Q, 1 )&0  TO  11  12400 

CALL  T 1 2L ' 1 .5)  12410 

I JKM-MODi I JK. 1 3 )  12420 

T I  ME 3  Y <  1  . JF )  ! 2430 

I F ( 1 JKM.GT .0 )G0  TO  10  12440 

WRITE  (6.700) IDEM  12450 

700  FORMAT  (1H112A&)  12460 

WRITE  (6.698!  12470 

10  WRI TE (6.657  IT IME. XYZ ( 1  ) . XVZ (2 1 .XYZ< 3 ) . WXS .TH, AH  12480 

WR!  TE(6.6961D(  1  )  .  0  !  2  )  •>  0  <  3  >  .  WV  S  .  V  12490 

WRITE (6.6  95  )B (  l  ) .8! 2 ) .8  <3  1  12500 

GO  TO  12  12510 

11  WRITE  (6.699)  12520 

12  IJK-IJK+)  12530 

RETURN  1 2540 

END  12550 

S1BFTC  TI2LP  LIST. REF. DECK  12560 

SUBROUTINE  T I 2L (KWIND. JF )  12570 

COMMON  V <22.5 ) «DY (22.5 1 .R0FF1  (22.5  i  12580 

C  12590 

C  IN  ALL  COMMON  STATEMENTS  D=DR 1 VER . As AOPTUN ,B=8FACT . E=EVAL . C  =  CHK . 1 = I NT  12600 

C  R2*RKERR2.R=RKG.T=TABL.F=PhAS1N.P  ° I NTR J « X=XOU T .AND  L=TI2LP  12610 

C  12620 

COMMON/XL/'  TH .  AH « D  <  3  1  •  8  (3  1  .  V  12630 

C  12640 

COMMONXOEL/  SLAT.CLAT.SLATG.CLATG.TlMEO.w  12650 

C  12660 

COMMON/DL/  RO  12670 

C  12680 

COMMON/ABXL/  XYZ (31  12690 

C  12700 

DIMENSION  A ( 3 . 3 ) . YD  <  6 1  12710 

C  12720 

C  IF  KWIND- 1.  THE  SECONO  DERIVATIVES  WILL  BE  COMPUTED  12730 

C  COMPUTATION  OF  THE  ROTATION  MATRIX  FROM  THE  EC  I  TO  THE  LAUNCHER  SYSTEM  12740 
C  W-R0TA1 ION  OF  EARTHSRS=RAUUjS  TO  LAUNCHER  SYSTEMSCLAT-COS ( GEOC  EAT)*  12750 

C  SLAT-S1NIGE0C  LAT )  CLATG-COS ( GEOD  LAT )  SLATG-S I N ( GEOD  LAT)  ARE  ALL  0  12760 

IN  DRIVER  12770 

12780 

WT«W» < V < l » JF )-T I MEO 1  12790 

A ( 1 . 1 1=  -COS< WT )  12800 

A ( 1 .2 ) =  ~SIN<  WT  >  12810 

A(1 .31=  0.0  12820 

A  { 2 . 3  1  =  CLATG  12830 

A ( 3  *  3 1 =  SLATG  12840 

A(2.1 )=-A(3.3)*A(l .2)  12850 

A  <  3 . 1  )  =  A(2.3t*A<1.2>  12860 

A(3.2)*-A<2.3)»A(1 .1  )  12870 

A (2 .2  )  =  A (3.3 )*A ( 1 . 1  )  12880 

12890 

COMPUTATION  OF  TRANSLATION  VECTOR  12900 

C  12910 

XL-RO*CLAT*A( 1 .2  1  12920 
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YL=-RO*CLAT*A  (  1  .  I  ) 

ZL*RO«SLAT 

C 

C  COMPUTATION  OF  POSITION  IN  LAUNCHER  SYSTEM 
r 

6(1  ) =  Y  <  1 7. JF )-XL 
S (2  >  =Y  < 1 8. JF ) -Yl 
B  <  3  >  = Y ( 1 9  «  JF ) — ZL 
CALL  MTRXL  (A.B.l.XYZ) 

C 

C  RELATIVE  MOTION  COMPUTATIONS 

c 

YD( 1  ) =DY ( 1 7. JF  >+W*Y < 1 8. JF ) 
vO ( 2 ) =DY ( 1 0, JF )-W*Y ( 1 7. JF ) 

YD(3)=DY( 19. JF) 

1 F ( KW 1 NO • NE • 1 )GC  TO  10 

YD (4  >=DY (20. JF  )+2.*W*YD (2 >+9*\«*Y < 1 7 . JF  > 

YD(5)=DY  <21 , JF )-2.*W*YD< 1 )+W*W*Y ( 1 8 . JF ) 

Y0(6)*DY(22.JF) 

C 

C  COMPUTATION  OF  VELOCITY  AND  ACCELERATION  IN  LAUNCHER  SYSTEM 
C 

CALL  MTRXL ( A. YD14 > . 1 .B ) 

10  CALL  MTRXL (A. YD( 1 ) . 1 .0) 

I F < KW I  NO « NE » 1  )G0  VO  20 

V=SORT  (0(1 )**2+0(2 )**2+0<3 )4*2 ) 

TH=57„2957795*ARSIN(D(3)/V) 

20  AH=57.2957795*ATAN2(D( 1  )  .0(2)  ) 

RETURN 

END 


sdata 

XNO 

1 

NASA  4. 

1  . 

AER08EE 

195  GAG  I 

TO-Z-LATi . 

4143.  4000. 

32.40392 

VEL 

283. 

0. 

1 

2.90 

348. 

ATMOS  404001, 

.00210161 100.99 

1 82 1 .09 

ATMOS 

450  1  , 

.00207931099,04 

1 787.51 

ATMOS 

5001  . 

.0020481 1097.09 

1 754.34 

ATMOS 

6002. 

.00198681093.18 

1689.76 

ATMOS 

7002. 

.00192681089.25 

1  627N03 

AT  MOS 

8003. 

.00186831085.31 

1 566.21 

ATMOS 

10005. 

.00175531077.39 

1 448,71 

ATMOS 

12007. 

.00164761069.4 

1341.01 

ATMOS 

1 4009. 

.00154  41061.36 

1238.69 

ATMOS 

160*2. 

.00144741053.25 

1142.77 

ATMOS 

19017. 

.00131  1040.97 

1010.27 

ATMOS 

22023. 

,001 18271028,55 

890.17 

ATMOS 

25030, 

.0010651 1015.96 

782.40 

ATMOS 

29040. 

.0009225998.96 

655 . 2 

ATMOS 

33052. 

.000795  981.65 

545.24 

ATMOS 
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